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INTRODUCTION 


With an assumed governing system of partial dif f erential equations 
(PDF's), numerical simulations typically consist of applying a grid generator 
and a PDE-solver. In each case, a communication link is established from the 
generated grid to the PDE-solver. This link occurs when the solver is written 
in terms of a grid which covers the physical region. In the event that a so- 
lu : on sh-uj'H develop large variations on a finer scale than is provided by 
the local g r ti spacing, a communication link must also be established from the 
r DE -solver 1 3 the g^id generator. The grid is then driven by the solution to 
meet the resolution requirements and is called adaptive. 

While the communication link from the grid to the solver is typically a 
matter of transferring solution data to the grid and then expressing the PDE's 
with respect to the grid, the link in the opposite direction is more subtle. 

To form it, we must use the solution data to determine the locations where the 
errors might become large. In a formal analytical sense, tne locations re- 
quK'trig resolution can be Detained by error estimates that are derived for the 
PPi or. I y^r- . Alternati vely f the same locations can often be determined direct- 
1 / from a basic knowledge of the problem undergoing simulation. From the di- 
rect viewpoint, we note that the solution is generally given as a vector cf 
dependent variables at each point in physical space and geometrically appears 
as a stir face over the physical space. Rather than taking the solution surface 
! t s ■■'•f, wo use our knowledge about it to extract the quantities which will ex 
perionce rapid variations. We then form another surface with those quantities 
and, to be descriptive, we call it a monitor surface since it must suitably 
monitor solution behavior. Rather than forming it with a direct vector of the 
quantities, scalar combinations are employed to achieve simplicity. The simp 
lest, situation occurs when a single scalar linear combination is used. By 
definition, the locations that require resolution are then determined by large 
gradients of this scalar. In a geometric spirit, the gradients are resolved 
when the monitor surface is uniformly covered with grid points. Moreover, the 
solution is more accurately represented by the grid when the bends in the 
morrtnr surface are also resolved. This can be achieved with curvature clns~ 
i er l ng . 

With the adaptive data suitably expressed in the form of an error esti~ 
mat. >•’ or a monitor surface, the next step is to develop grid generators that 
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can use the data to cluster points in the desired locations. The manner in 
which the clustering is performed must also be done automatically and reliably 
since human intervention is not assumed to be available for corrective ac- 
tions. In broad terms, the clustering is achieved by the movement of grid 
points and by changes in the number of grid points. 

Most methods which alter the number of points tend to start with a fixed 
global grid and proceed to locally add points as they are needed. Sometimes a 
means to remove points when they are not needed is also employed. As a conse- 
quence of local point additions and subtractions, the overall grid structure 
becomes rather complicated in the sense of the accompanying data structures 
and internal boundary treatment between fine and coarse spacing. In contrast, 
methods which move the points of a curvilinear grid maintain a simple data 
structure and have no such internal boundaries. However, with a fixed number 
of points, local resolution is achieved at the expense of depreciating the 
resolution in other regions. In many cases the depreciation is quite minor 
because the other regions often have an overabundance of points that would not 
significantly contribute to the accuracy of the simulation. The limitation of 
fixing the number of points then only occurs when there are not enough points 
to resolve either the local phenomena or the other regions. As a consequence, 
a reasonable strategy is to globally alter the number of points in each coor- 
dinate direction while maintaining grid point movement as the primary adaptive 
mechanism . 

One simple way to make the global decision is to first determine the 
number of points required for a mildly varying solution and then to estimate 
the number of extra points required for the severe variations. Along a coor- 
dinate curve, the estimate is just the number of severe variations times the 
number of points required for such variations. To adequately employ all 
curves in a given coordinate direction, the estimate to be used is then the 
maximum of the estimates for the curves. By adding it to the number required 
for a mildly varying solution, the total number of points in a given direction 
is obtained. With such strategies, the number of points can vary between the 
requirements of mild and arbitrarily severe solution behavior. The use of 
points is then somewhat optimal in the sense that there are just enough to re- 
solve everything to our satisfaction without a tremendous excess of modestly 
contributing points. 

To keep our discussion within simple bounds, we will only consider basic 


- 2 - 



grid point movement schemes with the implicit understanding that any of them 
can be applied as part of strategy where the number of points in each coordi- 
nate direction can also be adaptively adjusted as indicated above. As such, 
we will refer to the schemes as adaptive grids. This is consistent with the 
notion of grids as discrete coordinate transformations as opposed to meshes 
that can assume an unstructured format. 

The general topic of grid generation has been the subject of a number of 
conferences. These occurred at NASA Langley [1] in 1980, at Nashville [2] in 
1982, at Houston [»] in 1933, and at Landshut, West Germany [•*] in 1936. The 
latter is the start of a two-year international conference sequence. The 
topic also has appeared in a conference on adaptive methods in Maryland [5] in 
19 c 3 and has assumed a role of growing importance in the AIAA Computational 
Fluid Dynamics Conference sequence as well as in the international series on 
Numerical Methods in Fluid Dynamics, published in the numbered sequence of 
Lecture Notes in Physics by Spr inger-Verlag. 

In addition, there are general surveys by Thompson, Warsi and Mastin [6] 


.n i 982 , bv Thompson [7] in i98*J, by Eiseman [8] in 1985, and by Eiseman and 
F.rlebacher [9] in 1987. More specialized surveys have also appeared in the 
conferences at NASA Langley and at Nashville. Or, the topic of adaptive grid 
generation, the orevious reviews were given by Anderson [10] in 1983 and by 
Thompson [11] in 1935. The present review represents a fairly thorough cur- 
rent account of this active and important subject. 


jJEIGHT FUNCTIONS 


At. the root of adaptive grids is the desire to generate good grids with 
o‘->] is that contain equal amounts of a specified weight function. When the 
generation process is performed in physical space, the error estimates on 
gradients of the monitor surface are directly inserted into the weight func- 
tion along with other quantities such as monitor surface curvature or a back- 
ground distribution. When the process is performed on the monitor surface and 
the consequent grid is projected down onto the physical region, the gradients 
are automatically resolved by the use of surface arc lengths. With the weight 
functions for surfaces, the emphasis is on geometric quantities that allow u> 
to cluster points in a manner that Improves the discrete representation of the 
surface. These weights are applied against surface arc lengths, areas, and 
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volumes as opposed to the corresponding projected elements in physical space. 
As a consequence, the main clustering quantities in the weights are the vari- 
ous forms of curvature. This includes normal and mean curvatures for the 
basic geometry and geodesic curvature for boundaries. 

Regardless of where the grid generation process is performed, the basic 
formation and utilisation of weight functions remain the same: only the at- 

tracting quantities and the cellular elements are changed. In the most gen- 
eral form, the weight is a positive scalar function that incorporates quanti- 
ties for a background distribution, for arbitrarily attracting points, and for 
improving the structure and smoothness of the grid. The simplest and most 
commonly used form is a linear combination. The general linear expression is 
given by 


w=1+cM+c«M«+ + n y ( 1 ) 

w t m 1 c 2 n 2 c m .. m ) 

where the above quantities are represented with nonnegative functions M; that, 
have nonnegative coefficients c. to indicate the level of importance attached 
to them. We now think of the functions M,. as masses which then more strongly 
attract points when they are large. For the arbitrary attraction to a gr 3 di 
ent, it is the magnitude of the gradient or some normalized version of the 
magnitude. For a background distribution, it is a derivative of that distri- 
bution. For an improvement in the structure or smoothness of the grid, it is 
a pull of the grid towards orthogonality or towards conformal conditions. Al- 
together, these masses represent a collection of attributes that compete with 
each other and with the first term of unity which represents uniform condi- 
tions. The uniformity i3 actually achieved if all the coefficients should 
vanish. The uniformly distributed items are the elements on which the weight 
i3 applied. By making any of the coefficients arbitrarily large, the asso- 
ciated attribute can be made to be dominant and to thereby dwarf the effect of 
the uniformity term. In such a situation, the first term is viewed as only a 
guarantee that the weight function is never zero. In some circumstances, this 
term is dropped entirely. 

EQUI DISTRIBUTION IN ONE DIMENSION 

With the motivation to directly cluster points as the weight becomes 
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large, the Tirst consideration is to move the grid to get the cells to contain 
equal amounts or weight. This equally distributes the weight over the cells 
of a grid and is called an equidistribution process. Upon accomplishment, the 
equalization results in small cells to accommodate large weights and larger 
cells to accommodate smaller weights. Since the basic principles involved can 
be most simply displayed in one dimension and since most of the adaptive grid 
generation methods have their roots established in that context, we shall 
first consider one-dimensional developments and then proceed to the higher di- 
mensional versions. 

The Di ffe r ential Statement 

In one dimension, equidistr ibution occurs when the spacing between points 
is inversely proportional to the weight. Assuming a transformation from a 
curvilinear coordinate £ to the arc length s along a given curve, the rela- 
tionship is given by the differential statement 

wds = cdf; (2) 

where c is the proportionality constant. In the mapping, a uniform distribu- 
tion in ^ is transformed into a nonuniform distribution in s. The desired 
spacing is achieved because d£ is also viewed as a constant which then means 
that the spacing ds must shrink or expand to accommodate respectively large or 
small weights. 

The Local Integral Statement 

In the finite sense of grids, the appropriate conditions come from local 
integrations that are taken between the successive grid points. For an index 
i, the Integration is from £ to £ i + 1 and the corresponding s i to s iH . The 
consequent local integral statement is 



( 3 ) 



3 . 
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where 


r - £ . = constant 
u + 1 l 

The second part is the requirement for a uniform distribution in £ that now 
nas assumed the form of a g^ib with constant spacing. In an intuitive sense, 
the total weight over each cell is simply required to be a constant. 

The Direct Grid Statement 


While a variety of schemes could be constructed by evaluating the local 
integrals in various manners (i.e., distinct quadrature rules), the most corn - ' 
mon method is to assume a value of weight at the center and to use it as a 
constant over’ the interval. This approximation for the integral leads to the 
explicit grid statement 


w , (s - s.) - constant 
i + 1 i 


where typically 


(*0 


The effect of the weight on the grid points is now more directly evident than 
it was in the previous statements. Specifically, as w^j^is increased, s^ and 
s • 4 i approach each other. 


The Backward Global Integral Statement 


For a direct uncoupled relationship between grid point locations, global 
rather than local integrals are considered. When the weight w is given as a 
function of only the location s along the curve, a direct integration of Eq. 2 
leads to the transf ormation 
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^ ^min F(s)_ 

E ~ E . F(s 

max min max 


where 


( 6 ) 


s 

F(s) = f wdx 

s . 
min 

The derivation consists of first determining the proportional ity constant of 
Fq. 2 from the total integral as 

pfs ) 

c = - ~ (7) 

E - E • 

max min 

and then using the integral up to the current location s to get Ec. 6. 

Within the context of application, a curve has been given as a function 
of a parameter s and then the new parameter £ has also been expressed as a 
function of s. In the language of mappings, the functional relationships be- 
come the map from an interval of values s to the curve and the map from the 
same interval to the interval for E. The latter map, however, is clearly 
backward" with respect to the straight composition of maps which is from £ to 
s and then to the curve. As a consequence, the map of Eq. 6 will be called 
th® backward global integral statement. 

An implication from the backwardness is the need to invert the trans- 
formation in order to apply it. In terms of grids, the inversion is most 
simply accomplished by a local linear interpolation within a table of integral 
values for the previous grid points < a£ < .... < on the interval 

s . < 3 < 3 . With midpoint weights determined in the manner of Eq. 5, the 

min max 

integral of Eq. 6 is approximated by the trapezoidal quadrature rule 


F (a ) 

k 


k-1 

y w. 1 .(a. , 

. , l-e/p l-H 

i = 1 * 



( 8 ) 


for k = 2, 3 N and by 0 when k-1. When the previous points a are 

taken as the piecewise linear approximation to curve arc length, the spacing 
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increments are given by the Euclidian distance norm as 


a . 


i + 1 


P. , - P. 
i+1 1 1 


( 9 ) 


where P , P„ P,. are the previous grid points on the curve. Returning 

1 2 N 

to the transf ormat ion of Eq. 6, uniform spacing in E yields values of 
( j - 1 ) / ( N - 1 ) on the left-hand side and thus the interpolation problem 


' (S J ) 


N-1 F(a N } 


{ -c: 


for the internal points s- where j = 2, 3, .... N-1. The boundary points are 
already known and are and - a^- For each j, the value on the 

right-hand side of Eq. 10 is used in a search to find its position within the 
table of values from Eq. 8. The result is an index ci - m(j) for which 


F(a ) £ r(s ) < F(a ) Ml) 

m J 01+1 

with possible equality on the right if it is the last interval m =* N-1. Using 
the fractional amount of entrance into the interval 


Cl . 

J 


F(s.) ~ F(a ) 

J 5_ 

F ( a ) - F( a ) 
m+ 1 ra 


( 12 ) 


the new values of Sj 


are given by 


s . 
J 


a + 
m 




M3) 


and the corresponding new locations on the curve are given by 


Q. 

J 


P 

rn 


+ 


a . 
J 




(14) 


In the event that the transf ormat ion of Eq. 6 is analytically defined, 
the inversion can be more accurately executed by a variety of point iterative 
methods. Moreover, if the curve is also analytically defined, then the grid 
point locations can indeed be very accurately determined. However, in the 
general adaptive context, such analytical data is usually not available and 
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such high accuracy in the grid point locations is rarely critical. As a 
consequence, the most widely applicable inversion is the simple numerical ap- 
proach that was just presented here. 

The Forward Global Integral Statement 

As an alternative, the inversion problem can be avoided altogether by the 
construction of a forward transformation in place of the previous backward 
one. The basic assumption is that the weight w must be considered purely as a 
function of the new curvilinear variable E,. While the functional dependence 
can often be readily shifted by a change of variables, this viewpoint is im- 
portant because w is seen to vary with respect to the number of grid points 
rather than with respect to specific locations on the curve. In cases where w 
is explicitly defined as a function of E,, a change to s becomes a real compu- 
tational task that requires a global iteration to suitably converge. 

To derive the forward global integral statement, Eq. 2 is divided by w 
and is then integrated in the same manner as before. The result assumes the 
form 


s - s . G(C ) 

min _ 

s - s . ' = G(£ ) 

max min max 


where 


r 

GU) = \ 


E, . 
min 


w(c) 


and where the proportionality constant in Eq. 2 is given by 


S “3 

max min 

_ GU) 


(15) 


( 16 ) 


Using the midpoint values of Eq. 5, the integral in Eq. 15 is approximated as 
in Eq. 8 by the trapezoidal quadrature rule which here appears as 
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1 


(17) 


GU k ) - AC 


k-1 

I 


i»l 



By displaying the constant increment AC in front of the summation, the fact 
that C represents a counter of points Is accented. Since the simplest counter 
is achieved by taking the index i itself, it is common practice to assume that 
AC * 1 • When the forward transformation of Eq. 15 is evaluated at the uniform 
points c - i, the new locations are given by 


= s 


min 


max 


- s 


] C(J> 

rain'G(N) 


w here 


( 1 8 ) 


G(j) 


j-1 

l 


i-1 



f or j = 2 , 3, .... N. For j-1, the sum G(1) vanishes and the result is s, = 

s min ’ 

From the forward character of directly producing new parametric locations 
s each new grid point location Q i is directly obtained from another evalua- 
tion. If the curve is only defined by the previous grid points P. ] , P 2> .... 

P„ then the evaluation must be accomDlished by Interpolation. The simplest 
and most commonly used approach is to find the interval in the previous para- 
metric grid a, < a 2 < • • • < a N that contains Sj and then to use local linear 
interpolation. This proceeds by first searching to get the index m - m(j) 
such that 


with possible equality on the right if m « N-1. Using the fractional distance 
into the interval 


a . 
J 


3 . 

J 


a 

m+ 1 


- a 


m 


a 

m 


the new locations on the curve are given by 


( 20 ) 
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(21) 


Q 


J 


P ♦ 
m 



- P 

m 


) 


In contrast to the backward case, the local linear interpolation here occurs 
with respect to the previous parameter grid. 

Altogether, we could reasonably get the impression that the forward 
transformation is a much more efficient approach than the backward transforma- 
tion. This is certainly true when the weight function is explicitly defined 
as a function of the grid point counter £. However, in the adaptive context, 
the variations in the solution are present at distinct locations in the physi- 
cal region regardless of the grid employed for a simulation. As a conse- 
quence, the weights are really given as explicit functions of the spatial 
length s. This occurs because the weights are formed from solution proper- 
ties. In the application of the forward transformation, the weights are then 
only implicitly given as a function of £ in the form w(s(£)). The result is 
the requirement for an iteration. From a given weight function defined with 
respect to s, the evaluations w. over an initial grid of points s^ immediately 
defines a functional relationship between the grid point counter E, = i and the 
values w L . With this data, the forward transformation is then applied to pro- 
duce new grid point values s^. A return is made to the weight evaluations and 
the process is repeated. It is stopped upon reaching the limit of a conver- 
gence criterion. 


The Differential Equation Statement 


While the constant of proportional ity in the differential statement of 
Eq. 2 was explicitly evaluated in the global integral statements, it is re- 
moved by differentiation in the differential equation statement only to re- 
appear in the form of a second boundary condition. The first boundary condi- 
tion was the initial value for the global Integrations as applied to the 
simple first-order equations. 

For notational simplicity, a subscript by an independent variable will 
represent differentiation with respect to the variable. From Eq. 2, a 
E~dif ferentiation of the proportionality constant vanishes and results in the 
second-order differential equation 
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( 22 ) 


(WS^ - 0 


U po 


n expansion and division by w, it assumes the form 


3 + — s 

« * ^ 


0 


(23) 


where the weight term has been isolated. Accenting this isolation, the equa- 
tion can be written as 


; + as. = 0 

EE E 


(24) 


where 


a = 


E 


(25 


Conversely, when any equation of the form in Eq. 2H is assumed, the se- 
condary equation, Eq. 25, may be employed to get a weight function. The 

result is given by 


w = c^ e 


\ ad£ 


(25) 


where c, is the constant of integration. This weight is clearly equidistri- 
buted since all of the steps leading to the pair Eqs. 24-25 can be readily 
done in the opposite direction. The important consequence here is that Eq. 2» 
can be viewed as a basic equidistribution law. This form is called the dif- 
ferential equation statement since it is another statement which is equivalent 

to the others. 

An inverse form may also be given and can be determined either directly 
from the differential statement of Eq. 2. or from a formal interchange of vari- 
ables s and E in Eq. 23- In the case of Eq. 2, an. s-dif f erentiation of E</ w 
is seen to vanish. Upon simplification, the result is in the same form as in 
Eq. 24 but with a sign change and an interchange of variables. It also as- 
sumes the Poisson equation format for 1-D which is given by 


’S3 


W 


(27) 
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and which in higher dimensions has been extensively used for elliptic grid 
generation as discussed by Thompson [12] from a general perspective. 

The Direct Finite Difference Statement 


While various finite difference schemes could be employed to solve the 
previous statement in terms of a differential equation, a more immediate ap- 
proach is to use the direct grid statement of Eq. This is an approximation 

to the differential equation (Eq. 22) that is readily available. By writing 
Eq. : i for i-1 and i and noting that the same constant appears on the right- 
hand side, we have 


V / 2 


(s . 


i + 1 


- S. ) - W. j (S 

1 1 -Vo 1 


Vi* 


which can be rewritten as 


( 28 ) 


'iJ/ 2 S i*1 


(w. 




W . 1 ) s , 
l -V? ] 


W. j.S. , 
1-1 


* 0 


When combined with the boundary conditions 


(29 } 


S N ‘ s max 

we have a simple tridiagonal matrix to invert to get the interior values s^, 
s^, ..., The corresponding grid points along the given curve are 

then determined in the same manner as in the global forward Integral state- 
ment. In addition, the need for iteration is also required since the weights 
appear at unknown locations. Only the index i i3 known. This is discussed in 
the paragraph containing Eqs. 19~21. 

Mean-Value Relaxation Statement 


With a view towards direct extensions into higher dimensions, we are 
motivated to cast the tridiagonal finite difference -form of Eq. 29 into a 
statement of point relaxation that contains an intuitively plausible basis for 
such extensions. Without this view, there would be no compelling desire to 
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consider point relaxation since it is substantially more efficient to use the 
tridiagonal form in the sense of direct Gaussian elimination. In higher di- 
mensions, by contrast, the efficiency can be retrieved by multigrid methods 
[ 13 ] and by the use of parallel computing environments [ 1 * 0 . 

By solving for the center point in Eq. 29, a form that is immediately 
suitable for point relaxation arises and is given by 


s. 

1 






(31 ) 


In the Gauss-Seidel application, the values on the right-hand side use the 
available current information while the point Jacobi application ignores the 
current information in favor of the data from the previous sweep over the 
field. With either method of relaxation, however, the formula in Eq. 31 does 
not yield the desired intuitive basis since the evaluations of s are unaligned 
with the evaluations of w. By adding s i to each side of Eq. 31 and then di- 
viding the result by 2, the apparent intuitive obstacle can be overcome. Upon 
appropriately grouping terms, we arrive at 


where 


s . 
1 


w. 1 .S . 1 , 


V// 


V/2 


(32) 


s . 


3 . 
1 


^i +s i 


i 


and 


9 i-V 


Wi 


(33) 


In the present form, we can now intuitively look at the weights over the in- 
tervals on either side of s^. Assuming uniform weights over each side, the 
weighting centers are at the midpoints as given by Eq. 33 and are in corres- 
pondence with the weight values there. In a physical sense, the formula is 
readily interpreted as a center of mass formula for the interval from s^_^ to 

3 i + 1 * 

In the sense of pointwise movement, the center of mass form can be recast 
to explicitly display the movement. This is accomplished by subtracting the 



previous position 3. from each side of Eq. 32. The result is given by 


s.-s. « 
1 1 


w, .1 ,d. ^ w j d 




W. i. + w, 1 

1+V2 l J /2 


(3*0 


where 


d - 

+ 



n 0 


( 35 ) 


- s 


i J / 2 


are the maximum movement distances in positive and negative directions from 

S . • 

1 

The subsequent mapping to the curve is accomplished locally. The local- 
ness occurs because the move to a new cannot exceed one-half of the dis- 
tance from the old s i to its neighbors and s^ + 1 . If both neighbors are 

old locations, then the halfway restriction prevents an interchange of s^ with 
either s-^ or 3 ^ when their respective moves are also considered: at best 

they may simultaneously approach either midpoint of Eq. 33 . When the current 
values of s ;-1 and are immediately used upon availability, the halfway 

restriction is not required to prevent such an interchange. The important 
fact here is that with either method of point relaxation, the center point 
cannot move outside of the interval from to s i+1 at any stage of 

relaxation. This then means that the previous global search used in the 
direct Gaussian elimination approach would be replaced by a local search 
between the two adjacent intervals. Beyond this stage, the mapping to the 
curve proceeds with the same linear interpolation over the selected interval. 


The Leas t Squares Statement 


In contrast to the previous statements which evolved directly from the 
original differential statement, the evolution will now be considered to start 
from a global expression where the competition among the different locations 
appears simultaneously. The basic action comes from the minimization of the 
weighted sum of squares 
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( 36 ) 


N-1 

A * i-i W i + 1/2 ^ i + 1 ~ Si ^ 


The competition among locations corresponds to the distinct terms in the sum. 
By viewing the weighting coefficients as a sequence of positive constants, the 
minimum is determined by setting 


3A 

3s 


2w j _ 1/2 (s.- 




2w j+i/2^ s j+r s j^ 


0 


(37) 


and is thus seen to give the earlier tridiagonal system of Eqs. 29 - 30. 

When the weighting coefficients are considered to be evaluations of a 
single continuously defined weight function, the points of evaluation are seen 
to be only a function of the index. This means that the weight is only a 
function of the curvilinear variable £ which is then evaluated at uniformly 
spaced points. In practice, the points are typically chosen to be the index 
itself. 

As a matter of distinction, if the weight had solely been a function of 
the position s, then the evaluations at each index j would have varied with 
the current associated location s j . The result for the minimization equations 
would have been the inclusion of weight derivatives and a nonlinear dependence 
upon s. From an iterative point of view, the nonlinear effect would appear as 
a requirement to reevaluate the weights and/or their derivatives when each Sj 
assumes a new value. 


The Constrained Least Squares Statement 

As an alternative to the minimization with respect to the pointwise lo- 
cations 3 f , the previous sum of squares can also be minimized with respect to 
the spacing between points. In an intuitive sense, the contributions from 
each term will be separated from its neighbors and will more directly produce 
the direct grid statement of Eqs. 4-5. The added complication for this desir- 
able separation is the constraint that the sum of all of the spacing incre- 
ments equals the total length. 

With the spacing between successive points being denoted by 

h i+1/2 = 3 i+1 ~ 9 i 


- 16 - 



the previous sum of squares (Eq. 36) is expressed by 


N-1 

A = J-, w i+1/2 h i+1/2 
1 = 1 

and is subjected to the constraint 


( 39 ) 


N-1 


l h i*1 


i = 1 


= s - s . 
/2 max mm 


(40) 


For the minimization process, a Lagrange multiplier A is introduced to form 


N-1 


s - « - »n v , /2 - ( 


s - s . )] 
max min 


(41 ) 


from which we proceed to vary each h i +^/2 in the same manner as before. The 
minimum then occurs when 


3B 


3h. 


1 * 1/2 


2w i + 1 /2 h i+1 /2 


A = 0 


(42) 


which is readily recognized as the direct grid statement of Eq. 4 and which 
can be rewritten as 


h. 


i+1/2 2w 


i+1 /2 


(43) 


to express each increment as a function of the single unknown A. By substi- 
tuting the expression for each increment (Eq. 43) into the constraint (Eq. 40) 
of 3B/3A » 0, the Lagrange multiplier is found to be 


A 


2 ( 


3 

max 

N-1 


l 

i=i 


w 



J 

i+1/2 


(44) 
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Returning to Eq. 43, each increment is now known. Upon also inserting Eq. 38 
into Eq. *43 , the k-th increment is then given by 


s 


k+1 


1 

W k»1 /2 (- 

N-1 , S max 


u w 

i-1 i + 1/2 



( 45 ) 


5y forming partial sums starting from k = 1, the forward global integral 
statement of Eq. 18 is now readily apparent. 

The solution to the constrained least squares problem also has a direct 
geometric interpretation that was observed by Steinberg and Roache [153. With 
the increments considered as the Cartesian coordinates of (N-1 )-diraensional 
Euclidean space, the basic sum of squares in Eq. 39 represents a hyper-ellip- 
soid for each value of A. Clearly, each such Cartesian coordinate h i+1/2 is 
bounded by 


l h i + 1 /2 1 " vw . 


i + 1 / 2 


(46) 


where equality is achieved only along the h^ + ^/ 2 ax ^* 

From this we see that the family of hyper-ellipsoids varies from the 

origin when the "radius" /a is zero to arbitrarily large sizes as the radius 

is increased. Above the origin, however, we have a hyper-plane of constraints 

defined by Eq. 40. Geometrically, it passes through each axis at values of 

q - a . Since hvDer-ell ipsoids are curved convex surfaces, there then is 
max min* 

only one hyper-ellipsoid that intersects the hyper-plane at a single point. 
This also is a point of tangency. While lesser values of radius /k have no 
intersection, higher values have infinitely many. But to minimize A, we are 
only interested in the smallest value with an intersection. This is just the 
single point where tangency occurs. 

Variational Statements 

The natural framework to continue the previous least squares statements 
is provided by the variational statements. In distinction, these formulations 
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are now considered from a continuum rather than a discrete point of view. The 
general statement here is to find the extremum of an integral. When the inte- 
gration is with respect to the curvilinear variable E, we shall consider the 
integrand type 

F 1 - [w(E)sJ a (*^7) 

for a weight function dependence on E and the integrand type 

F 2 > [w(s)s^] a (W 

for a weight function dependence on s. When a - 1 for F 1 in Eq. *17, the inte- 
gral becomes the continuum version of the least squares statement. Upon ap- 
proximation with unit increments in £ , the relationship becomes 


N N-1 

I - 1 Vl/3 (s i.1 ■ 3 i> ! 

1 ^ 1=1 


(1*9) 


and provides a direct connection with the previous discrete least squares de- 
velopments . 

Returning to the integral forms j F^E, the extremum is found as before 
by setting first derivatives to zero. This is accomplished by using the Cal- 
culus of Variations [16]. The result is the Euler equation 


( 


3 _ 

3s 


dE 3s 


: )F i 


(50) 


By direct computation, it can be readily verified that the Euler equations re- 
duce to a differential equidistribution statement for each choice of inte- 
grand. In the case when the weight function depends upon the location s, we 
also note the appearance of a weight derivative. This was also observed when 
we examined the basic least squares statement. 

As an alternative, we may change the differential in the integral by us- 
ing dE = E ds. We mu3t the 11 , however, also use a = i in Fj. By the same 
argument applied to the interchanged variables, the Euler equation is now 
given by 
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(51) 



d_ _a_ 

1 = 


>< F l« a ) 


0 


By carrying out the algebraic manipulations for each given F^, the differen- 
tial equidistribution statement can be retrieved as earlier. The selection of 
integration variable is then clearly seen to be a matter of convenience. 


Evolutionary Statements 

In the previous statements, equidistribution was expressed in various 
spatial forms that produced a distribution either directly or by iteration. 

The need for iteration arose in situations where the weight, being known as a 
function of s, was instead used as a function of I * i rather than s. This 
occurred in a natural fashion when the weights were needed at yet undetermined 
positions s^ the only possibility was then to use the currently available 
values in an iterative cycle. During the course of iteration, no considera- 
tion was given to the potential control of the evolutionary iterative cycles. 
In the present section, such an evolutionary control is examined. 

In the most immediate manner, the evolutionary control can be stated by 
the parabolic partial differential equation 

Ks - (ws) (52) 

U >■» s 

where the non-negative function K = K(£) represents the rate of evoluation for 
each point £ = i. This rate varies from an infinite speed when K vanishes to 
a zero speed when K is infinite. As a consequence, the motion can be greatly 
reduced or stopped in certain locations and can be gradually increased to 
large values in other locations. From a physical point of view, the diffusion 
of points into the equidistribution positions increases with 1/K. 

In the case where K vanishes everywhere, the points then instantly move 
into the locations for the equidistribution of w. The case corresponds to the 
earlier differential equation statement of Eq. 22 and its various forms in 
Eqs. 23-27. The infinite speed represented by the instantaneous adjustment 
comes from the elliptic nature of the consequent differential equation. In 
practice, however, the elliptic type of differential equation would still have 
to be solved numerically and, in the given situation, would have to be done by 
iteration. The instantaneous adjustment is then represented by the direct mo- 
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tion into the converged positions. In the iterative cycle there is simply a 
succession of finite speeds. This can be seen, for example, by the mean value 
relaxation form in Eq. 3^- 

With a spatial grid, the parabolic partial differential equation of Eq. 

52 assumes the form 


ds . 


K. 


i dt ’ W i + 1 /2 3 i+ 1 



w 

1-1 / 2 




(53) 


where the grid point velocity at i is explicitly displayed. In an algorithm, 
the positive numbers w. +1/2 and w jL _ 1 /2 are first computed and then the 3j_ is 
moved accordingly. When the original s^ lies between s i _ 1 and s i+1 , both the 


•1/2 term w. +i / 2 (s. + 1 ~'Sj) and the i-1/2 term w i _ 1 / 2 ^ s i 


s._ 1 ) are posi- 


tive. If the 1*1/2 term exceeds the i-1/2 term, then s^ must move in the 
positive direction towards s i+1 to simultaneously shrink (s i+1 - Sj) and ex- 
pand ( s - 3 ). If the i * 1 / 2 terra is less then the i-1/2 term, then the 

velocity is in the negative direction towards s i _ 1 . If both terms are equal, 
th°n the forces balance and thus there is no motion. When the original s^ for 
an algorithmic step lies outside of the interval s < s < s J+1 , the motion 
is directed back towards the interval. For example, if the original s^ ex- 
ceeds s,- + 1 , then both terms are negative and thus the velocity of eq. 53 is 
also negative, sending s^ towards or past s^ +1 - Moreover, when the ordering 
of the points s^ is monotoni cally decreasing rather than increasing as in the 
discussion just given, the same conclusions result with only a change in 
directions . 

As a further observation of Eq. 53. it is noted that when the velocity is 
zero, the tridiagonal system of Eq. 29 is obtained. That system is rapidly 
approached when the diffusion rate 1/Kj is sufficiently large. In the evolu- 
tion towards the* earlier tridiagonal system, a temporal discretization must be 
considered. The simplest is the single step from time t n to time t n+1 to- 
gether with a backward temporal difference over the time step At - t n+1 - t n . 
With the time indicated by a superscript, the parabolic partial differential 
equation is approximated by the fully implicit system 


n+1 


,,n f i M - u n fs n+1 - * n+1 ) - „ n U n 

V At J i+1/2 l i+1 i J i-1/2 1 i 


n+1 n+1 i 

" 3 i-1 ] 


( 54 ) 


which assumes the tridiagonal form 
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( 55 ) 


n n+1 
w i + 1 /2 3 


n+1 e n 
i + 1 ^ W l + 1 


/2 


w + 

i-1/2 At 


-Is" 

t j3 i 


n+1 ♦ W n 3 n+1 

1 - 1 / 21-1 


K? 

— s" 
At 3 1 


Upon Inspection, strictly positive values of k" ensure diagonal dominance. 
Moreover, very large values will force an essentially diagonal form and there- 
by cause the i-th point to essentially be frozen: that is, an essential re- 

At the opposite extreme, a zero value of will give 


... . n+1 

duction to s. 

1 


i 


the earlier tridiagonal form of Eq. 29. As from the earlier form, the same 
manipulations will lead to a mean value relaxation statement. With a point- 
wise update from only n-level information (point Jacobi relaxation), the 
statement becomes 


n+1 /2 
3 i+1 /2 


K n 

n n i n n n 

W i + 1/2 S i + 1/2 At S i Vl/2 3 l-1/2 

K? 
i n 

W + <5 + u 

i+1/2 At i i-1/2 


(56) 


and reduces to Eq. 32 when K™ vanishes. In continuation, the least squares 
statement of Eq. 36 extends into the form 


A 


N-1 


l 

i-1 


n+ 1 


{KjAtf- 


At 


+ 


n r n+1 
w i+1 /2 S i+1 



(57) 


and reproduces the tridiagonal form when it is minimized with respect to var- 
iations in s"* 1 . Here, of course, K 1 = 0. The pattern of development here 
clearly proceeds through the more general variational format. 

Further accuracy in representing the parabolic partial differential 
equation can also be developed in a similar manner. This would come from 
writing the finite difference equation about the (n+1/2)-time level in a 
Crank-Nicolsen form. The result would be second order rather than first order 
in time provided that the weights are appropriately linearized about the time 
level t n . The associated tridiagonal form would then have more complication 
in the generation of its coefficients and source term. As a consequence, the 
simpler first-order accurate form is the likely preferred form. 
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THE SPECIFICATION OF COEFFICIENTS IN THE LINEAR HEIGHT 


Fractions 


When linear weights are employed for equidistribution, there is a dis- 
tinct advantage in using the backward global integral statement as opposed to 
some of the others. The advantage is the capability to more precisely specify 
the level of importance of the various clustering quantities in the weight. 
This comes from the linearity in both the integration and the weight. With 
the general linear weight of Eq. 1 , the integral of Eq. 6 becomes 


m 

F(s) = H q (s) + l c.H.(s) 
i = 1 


( 58 ) 


where 


Hq(s) s s min 


and 


H i< s > • I Vx 

min 

Upon evaluation at s = s max , the equation becomes a relationship between total 
amounts. With the total length 


L = H 0^ s max^ = 3 raax 3 min 


( 59 ) 


and the total amount of each quantity 


the total weight integral is given by 


( 60 ) 


F(s ) - 
max 


L + 


Vi 


c I 
m m 


(61 ) 
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Since each term is positive and represents a part of the total weight, a divi- 
sion by the total weight results in the fractional decomposition 


f 0 * f 1 


‘ m 


( 62 ) 


where 


‘0 F ( 5 ) 

max 


is the fractional contribution from the total length and 


(63) 


f . 

i 


c.I. 

l i 


F ( s_ 


max 


(6H) 


is the fractional contribution from the total amount of the i-th quantity. 
Since each quantity must have a specified level of importance c^ and since a 
direct specification would be somewhat ad hoc, a more precise strategy is to 
specify the fractions of Eq. 6K and then solve for the consequent c^ With 
these specifications, the decomposition of Eq. 62 produces a value for the 
fractional amount of length f Q which then can be inserted into Eq. 63 to yield 
the total weight integral in the form 


F(s 


max 



(65) 


Upon substitution into Eq. 6**, the weighting coefficients are then given in 
terms of specified values as 


c. 

l 




f Jl. 

m l 


( 66 ) 


When these coefficients are employed in the integral F(s) of Eq. 58, the cor- 
responding backward global integral statement of Eq. 6 becomes 


£ - 
^max 



m H. (s) 

y r — 

A i H.(S ) 
i=0 i max 


(67) 


where f Q has now been reinstated. Upon observing the form of Eq. 67 together 
with the argument leading up to it, the role of fO is seen to be assumable by 
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any of the other fractions appearing in Eq. 62: that latter expression merely 

provides a partition of unity which determines any fraction when the other m 
fractions are specified. The fraction so determined then plays the same role 
as Tq above. To guard against the possible indeterminate situation should 
Hi ( 3 max> vanish * a 3afe practice is to replace each H 1 (s max ) by H j[( 3 max ) + e 
for some very small positive number e. The development of the above fraction- 
al specifications was originally derived by Dwyer [17] for cases up to m * 2 
fractions and was then extended by Eiseman [8] to cover cases with any number 
of fractions, m. 

Level of Significance 


When the consideration is shifted from the proportion of the weight 
F(s ) to the proportion of the length L, the specified fractions above must 
undergo an adjustment to account for multiple effects that occur at the same 
location. The spirit of the adjustment can be readily examined in the case of 
just one effect since it must then compete with the uniformity term that pro- 
duced L. With m = 1 , the first step is to assign a level above which Mj is 

considered to be significantly large and thus important enough to be counted. 

The effect here is to remove the minor level activity in M 1 that is typically 
spread over large regions and to thereby emphasize only the local regions of 
interest. The assignment of a then defines the set of local regions as 

(s: M 1 > } . The Lebesque measure [18] of the set represents the length 

over which is significant. With 

« meas{3: } (68) 

the length where is unimportant is up to the order just (L-L^). By con- 
trast, the total length where M 1 is important is counted twice: once by L 1 

and once by c^I^ in correspondence with 1 and in w. As a consequence, 

the total weight is partitioned in the form 

• a -V * ( c ,V L i> (69> 

and fractions are now defined by 
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( 70 ) 


F(a ) 
max 


and 


vrs 

F(s max } 


(71 ) 


With f Q - 1-f 1? the equations are solved as before to produce the weighting 
coefficient which now becomes 


f.L-L. 




(72) 


The need for the protective e in Eq. 67 can be dropped here since a vanishing 
1^ would mean that L| * L and hence it would then be reasonable to set c^ * 0 
thereby removing the effect of M^ altogether. On reflection of the above pro- 
cess, one may also note a parallel resemblance to studies in the area of sta- 
tistics where it is common practice to assume a level of significance [19] for 
events . 


Minimum and Maximum Spacing 

When there is only one coefficient to be specified, its specification can 
be translated into the more intuitive and geometric specification of a rela- 
tionship between the actual minimum and maximum spacing. This is derived **~om 
the direct grid statement of Eq. 4 by expressing it under both minimum and 
maximum conditions. The result is given by 


w (As; . - w . (As; 

max min min max 


'73) 


When the general linear weight of Eq. 1 is written for a single coefficient (m 
^ 1 ) in the form 


w » 1 + cM 


1 Tl| ) 


and m inserted into the relationship of Eq. 73* that coefficient is then de- 
termined by 
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c - 


R - 1 


M - M . R 
max min 


where (75) 


R 


(As) 


max 


(As) 


min 


As a consequence, the ratio R of the actual maximum to minimum spacing can be 
specified and will determine a valid coefficient c provided that 


R < ^ (76) 

mi n 

This condition merely states that the global variation in specified spacing 
cannot exceed the corresponding variation in the weight. The specification of 
such ratios R is due to Nakahashi and Diewert [20]. Upon observing the above 
development, it is also worth noting that the same argument can be applied to 
the general linear weight when only one coefficient remains to be determined. 


THE ATTRACTION TO A GIVEN GRID ON A CURVE 


In a number of circumstances there is a need to have some adherence to a 
specified distribution of points along a curve. These needs will typically 
arise, for example, when an expanding grid is desired for far-field condi- 
tions, when a nearly orthogonal grid is desired for enhanced grid structure, 
and when a smooth temporal evolution is also desired. While cases such as the 
expanding grid can be given by either specifying a weight function or a grid, 
the structural and temporal properties appear only in the form of a grid. The 
usual grid for the structural aspects comes from the computation of orthogonal 
trajectories [21]. They evolve from the adjacent and most current coordinate 
curves to our given curve which here is assumed to be a part of a higher di- 
mensional grid. The grid for the smoothness in temporal evolution comes from 
the previous grid or grids in a time-like sequence of grids. 

Rather than considering further situations in which it is desirable to 
adhere to given grids, we shall next consider how to accomplish such adher- 
ence. This will proceed first in a direct fashion by inverse equidistribu- 
tion, then by diagonal dominance, and finally in several variational forms. 
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Inverse Equidistribution 


The most direct method of providing attraction to a given grid i3 to per- 
form the equidistribution process in reverse. That is, we derive the weight 
function from the given grid. The result is an approximation to the deriva- 
tive of the desired distribution function had it been available in a smooth 
continuum form rather than as a grid. Upon applying the weight in a direct 
sense, the specified grid will be reproduced up to errors in the derivative 
approximation and in the subsequent quadrature for the direct equidistribu- 
tion. Assuming that such errors are sufficiently small, the weight derived 
from the grid is then viewed as only an attracting mass to be inserted into 
another weight function where it must then compete with other attracting 
masses. In terms of the general linear weight of Eq. 1, the mass is one of 
the functions M^. As the corr esponding coefficient Cj_ increases, the result- 
ing grid approaches the approximate reproduction of the specified grid. 

The forces in a weight function, however, are applied to directly control 
the spacing between the grid points rather than the actual positions of pai 
ticular grid points. In the general situation, it is the local concentrations 
of grid points which is achieved, not their actual positions. For any local 
concentration, grid points are attracted from the entire curve to satisfy the 
local spacing requirements. This local attraction will cause a shift to ap- 
pear in all points. As a consequence, that shift will directly appear in the 
situation when the local attraction is balanced against an attraction to a 
prescribed grid, tfhile the spacing in the prescribed grid will be rendered in 
a fair and competetive manner, the shift will lead to unacceptable results for 
the actual positions. The problem here is that the attraction to a specified 
grid is position-sensitive while the equidistribution process is position-in- 
sensitive. The primary cause of the problem is that local constants of inte 
gration are being effectively introduced from each local cluster. The only 
apparent remedy is then to employ an iterative scheme about the equidistribu- 
tion process. At each stage, an equidistribution of the current weight would 
produce results that in turn would lead to a new weight for the next cycle. 

In this framework, the inverse equidistribution strategy provides a distribu 
tion function that would enter the iterative procedure. Although such a pro- 
cedure has not yet been developed, the process of inverse equidistribution 

has . 
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To more closely examine the inverse equidistribution process, we shall 
explicitly derive weights at the current grid points 


s 1 < s 2 < . • • < s N (77) 

for the purpose of providing pointwise movement towards a specified grid 

s 1 = r l < r 2 < • * ’ < r N = S N (78) 

From a continuum viewpoint, the inverse equidistribution is accomplished by 
setting 


w - (79) 

for then the constancy of ws becomes the constancy of s p and hence a linear 
relationship between s and r that becomes an equality since s 1 = r 1 and s N = 

r N * 

From a discrete constructive viewpoint, we invert the direct grid state- 
ment of Eq. to get 


1 

w. « ■ 

1,1/2 r i*1 - r l 


(80) 


from the specified grid of Eq. 78. The unity in the numerator was chosen for 
convenience: any other nonzero constant would produce the same result since 

it would disappear in the ratio of integrals in Eq. 6. With Eq. 80, the 
weights are established at r 1+1/2 and must either be used there or be trans- 
ferred to the current points s^. The latter case may proceed either directly 
or with an intermediate transfer to points r^. In the direct approach, the 
derivative profile is obtained by an interpolation between half points n i+1 / 2 
and an extrapolation to end points r 1 and rjj. 

The simplest profile is given by a piecewise-linear interpolation that is 
extended by a horizontal extrapolation attached to either end. The piecewise 
linear transfer of weights is accomplished by searching the points 


r 1 < r 1/2 < < r N-i /2 < r N 


(81) 
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to find the interval that contains s k for each k. If the interval Is either 
the leftmost or rightmost one, then the weight is respectively either w 1 /2 or 
w N-l/2 as determlned Ed * 8o - otherwise, we have r m _ 1/2 S s k < r m+1/2 for 
some m = m(k) and hence the weight 


w + — 

m-1/2 r 


3 k r m-1/2 


ra + 1/2 r m - 1/2 


( W m+1/2 ‘ W m-1/2^ 


( 82 ) 


This weight can then be used as an attractive mass at the point s k within the 
linear weight of Eq. 1 . 

To employ it effectively, the values at s k must represent a fair sampling 
cr else the basic features of the weight could be lost or distorted. A poor 
representation here is certainly likely to increase the quadrature error upon 
application. This circumstance can readily arise if there are too many points 
s in any given interval from Eq. 81 . There would then be regions where the 
weight is represented too linearly. As a consequence, it Is advisable to use 
weights from Eq. 82 only when the current (Eq. 77) and the specified (Eq. 78) 
grids are nearly interlaced. Even still the inaccuracy may be too large. 

This would then suggest an application only when the r^ and s^ are nearly 
equal . 

Being somewhat forced to consider cases where the s^ and r^ are not too 
far apart, a further option has arisen. When the backward global integral 
statement of Eq. 6 is being employed, the curve together with the distribution 
function is being represented at existing locations Sj of Eq. 77. Upon in- 
spection, those locations could readily be replaced by the locations associ- 
ated with the specified grid for r^ given by Eq. 78. With that replacement, 
the quadrature of Eq. 8 produces a backward global transformation from Eq. 6 
which precisely moves s^ into r^. Linder such a replacement, other weighting 
quantities together with the curve description must undergo a transfer to the 
new positions r^ and thus represents another source for some error. This 
leads to the enrichment strategy whereby the points rj are merely included 
along with the original Sj^ in the curve representation. In essentially the 
same manner, the exact reproduction of the given grid locations r i is pos- 
sible. The derivation represents a modest extension. These developments 
leading to weights that reproduce specified grids were jointly constructed 
with this author and Michael Bockalie. 
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Diagonal Dominance 


In contrast to the inverse equidistribution approach above, we may place 
the attracting mechanism outside of the weight rather than in it. With this 
objective in mind, we return to the direct finite difference statement of Eq. 
29 and then add D^Sj-r^) to the right-hand side. With the assumption that 
the numbers Dj are nonnegative, we obtain the diagonally dominant tridiagonal 
system 


“i+l/2 3 i+l 


- (u 


1 + 1/2 


Vl/2 )3 1 


W i-1/2 S i-1 


- D.r i 


( 83 ) 


which as before is fully defined when boundary conditions such as those given 
by Eq. 30 are assumed. As increases, the diagonal dominance is intensi- 
fied. In the limit, the system approaches a simple diagonal form and this 
causes each current grid point s i of Eq. 77 to approach the corresponding spe- 
cified grid point r^ of Eq. 78. 

By providing an index to D^, we have included the possibility to approach 
the specified grid at a nonuniform rate: a uniform rate would come from set- 

ting D^ = D. As a consequence, we may also selectively approach only parts of 
the specified grid by smoothly adjusting the values of Dj from zero to local 
maxiraums at the desired locations. The smoothness is viewed by considering a 
function D of E to be smooth: the values at £ - i being denoted by D^. 

Moreover, we may also consider a sequence of specified grids together 
with a corresponding sequence of such functions. By adding the associated se- 
quence of terms to the right-hand side of Eq. 29 and following the above rea- 
soning for a single term, the simultaneous attraction to the sequence of grids 
can be considered. The relative sizes for D(E) in each term give the relative 
importance of that term. In the application, there i3 then a balancing of im- 
portance for the attraction to the separate grids. This is the same type of 
balancing operation as in the previous case with inverse equidistribution 
where each specified grid was represented by an. attracting mas3 within the 
weight. 

To explictly examine the relationship between the equidistribution and 
diagonal dominance approaches, we consider the case of attraction to one pre- 
scribed grid of points r A as displayed in Eq. 78. In the diagonal dominant 
approach, the attraction is represented by Eq. 83+ Since the weights are 
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functions of s rather than i and since the positions represented by s^ are 
unknown, iteration is required. As in the earlier tridiagonal approach, 
weights are evaluated from current values of position s^ and then the tridi- 
agonal system is inverted to obtain new locations Sj from which the process is 
repeated. Now let the fixed location r^ in Eq. 83 be given by the relative 
location 


r . 
1 


a i 3 i+1 /2 


(1 


a i >3 i-l /2 


( 8*0 


through a^. The half-point values are the same as in Eq. 33 and are deter- 
mined from the current locations in the iterative cycle. By inserting the 
local repr esentation of Eq. 8*4 into the diagonally dominant tridiagonal form 
of Eq. 83 and regrouping terms we get the new adjusted weights 


w 

w i + 1 12 
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i-1 / 2 


w + 

i + 1 /2 


w . 


i-1 /2 



. D . 

1 1 

- a . )D . 
1 1 


(55) 


for the original tridiagonal form of straight equidistribution as given in Eq 
29. While this yields a local equidistribution statement for each i, the 
global form is somewhat altered. That is because the half-point weight eval- 
uations from Eq. 85 will vary with the equation. In particular at i+1/2, the 
application of Eq. 85 yields 


w , . - w . „ + -r a . D , 

i + 1/2 1 + 1/2 2 1 1 


and 


( 86 ) 


w = u +• 

i+1/2 1*1/2 


r(1 “ 


a i-H )D i+1 


respectively for equations about i an d 1+1. To employ the weights in any of 
the quadrature forms presented, we need uniqueness and this leads to the re- 
quirement that 


i +1 




(87) 


which then determines the diagonal scaling up to a single multiplicative con 
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atant . With this determination, we then immediately observe that any of the 
basic forms for equidistribution can be applied within such an Iterative cycle 
of global updates. Because of the specialized choice for Dj in Eq. 87, it is 
clear that the diagonal dominance approach has some added and useful general- 
ity. The freedom to choose the form of D means that the attraction to the r A 
can be locally controlled. 

The Variational Format 


In correspondence with the previous forms of attraction to a prescribed 
grid, there are variational statements both with the attraction in the weight 
and with it outside of the weight. With it in the weight, we can simply use 
reciprocal derivative of Eq. 79 in the integral of Eq. 49 that represents the 
continuum version of the discrete least squares statement. This was the 
starting point for Steinberg and Roache [15] who considered such specified 
grids to be "reference grids" from which various metrical properties could be 
extracted and then directly used to form weights. 

With the attraction outside of the weight, the minimization of 

j [ws z +■ D(s-r) z ]d£ (88) 

provides the same form as in Eq. 83 with the same function D(£) that controls 
the attraction of s to a specified r. The r assumes the previous role of a 
given grid and is thus a function of the spatial location s rather than the 
eventual grid point counter £. This clearly penalizes the basic equidistribu- 
tion in favor of attraction to r as D increases in size. In a similar spirit, 
this statement can also be applied to derivatives. In particular, we have the 
form 

j [ ws z + D(s^-r^) z ]d£ (89) 

that was considered by Bell and Shubin [22] to enforce temporal smoothness. 
This was achieved by taking r to be the previous time level of s in the solu- 
tion of a one-dimensional time-dependent PDE. In the application, the weight 
was also smoothed: this being done by a spatial filter . With D as a con- 

stant, the Euler Equation reduced to the form 
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ws„ ♦ D(s-r) r * constant 


(90) 


which then can be simply integrated. This can be contrasted with the more 
complex integration that comes from Eq. 88. 

EVOLUTIONARY FORCES 

As in the discussion on the various forms for equidistribution in one 
dimension, we are now in a position to examine some of the forces which can be 
applied to some advantage in the evolutionary setting. The forces will appear 
both in the context of equidistributed weights and outside of such weights. 
Within the weights, the local application from a global determination will be 
witnessed. Outside of the weights, the attraction to a specified grid will be 
examined, the minimum and maximum spacing controls will be displayed, and the 
control of cell eccentricity will be explored. The role of diagonal domi- 
nance, here, is seen to arise again in a combined sense. 

Globally Determined Weights 

In the earlier consideration of weight functions, the evaluations at 
half-point locations were typically assumed to come from only the nearest 
neighbors as given in Eq. 5 of the direct grid statement for equidistribu- 
tion. Rather than using only that simple local average, we shall consider a 
global determination for the local weights about the point £ - i for applica- 
tion about the corresponding s^. For this purpose, suppose that we have es 
tablished some error indicator at each point. By using an average error 
value such as 
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the half-point weights about i can be written in the form 
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where c is a prescribed positive scaling constant and a Is an attenuation con- 
trol. The effect of a is to give more importance to the error deviations 
which are closer to i. As a increases, the influence of the local error grows 
while the more distant points are more attenuated. When the weights of Eq. 92 
are inserted into the evolutionary equidistribution statement of Eq. 53. we 
obtain the grid velocity equation. 
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This equation represents a local equidistribution since the weights about i-1 
and i 1 will produce generally distinct values at the respective half points 
i-1/2 and i -*- 1 /2 . The circumstance is just a parallel to the case with the 
diagonally dominant attraction to a given grid as represented in Eq. 86. 

With the evolutionary equidistribution of Eq. 93. we are very close to 
the grid velocity model that was proposed by Rai and Anderson [23] on the 
basis of physical intuition. They formed their model on the basis of an 
analogy with gravitational forces. That, for example, is where the attenua- 
tion from an inverse power law in distance from i was conceived. In dis- 
tinction from what is presented here, their model is obtained from Eq. 93 by 
replacing both (s i+1 " ) and ~ ) with a single factor such as 

(a - s^_^ ) and by combining K^/c into a single prescribed constant to 
uniformly scale the right-hand side. 

Returning to the choice of weights in Eq. 92, it should be noted that the 
essential character is the attenuation with respect to distance from the loca- 
tion of application. It is not the particular choice of functional form: 
many other forms could be employed for the same purpose. The same comment 
also applies to the implementation of the error indicator: the basic feature 

is that it can change sign. While thi3 would be undesirable in the nonevolu- 
tionary context, it is beneficial here and provides and additional restoring 
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force. 


The Attraction to a Given Grid 


With the understanding that a given grid such as that displayed in Eq. 78 
comes from a given distribution function r(s), the basic evolutionary para- 
bolic partial differential equation of Eq. 52 can be appropriately adjusted to 
provide an attraction to the fixed r(s) and to the fixed grid of r 1 values. 
Under the adjustment, it assumes the form 


Ks = (ws) + D(s-r) 


(9H) 


where the attraction is represented by the scaled quantity s-r. When the suc- 
cessive steps leading up to the fully discrete form in Eq. 55 are repeated, 
another tridiagonal form is obtained and is given by 

n n+ 1 r n n i n + l n n+1 

V1/2V1 " ( Vl/2 + Vl/2 + It + °iK + W i-1/2 S i-l - 
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In a combined sense, diagonal dominance is achieved from both the temporal as- 
pect and the attraction to r^. The respective effects, however, are essen- 
tially separated: as s L converges, the influence of decreases while that 

from Dj remains. The reason is because the terms converge towards an exact 
cancellation while the terms generally do not. On comparison with the 
earlier diagonally dominant attraction of Eq. 83 , the evolutionary form pre- 
sented in Eq. 95 represents a natural extension. 


Minimum and Maximum Spacing 

At this stage, the use of minimum and maximum spacing has been developed 
only for a linear weight with only one constant to be determined. That deter- 
mination appeared only in the form of a ratio between minimum and maximum 
spacing and was given by Eq. 75. Neither the actual spacing nor the multiple 
constant environment were considered. In contrast, such a consideration can 
be readily developed in the temporal setting. 


- 36 - 



To develop the required forces, the definition of the parabolic partial 
differential equation is again extended by the inclusion of new terms on the 
right-hand side. The generic extended form is given by 


Ks - (ws + A) + B 

L S 


( 96 ) 


where forces over the velocity can be inserted as terms in either A or B. In 
general, A and B represent sums of the forces which do not readily fit into 
the equidistribution of w with respect to s. While the position of A appears 

to be somewhat close to that equidistribution, it gives us the opportunity to 

rationally create similar forces without considering the factor of some As as 
would appear in the weight w. The location B allows for an even more direct 
insertion of force as has already occurred in Eq. 9^ when the (s~r) term ap- 
peared. That term provided an attraction of s^ to r^. In the overall cir- 
cumstance, the various forces in w, A, and B all act in a concerted fashion to 
provide a grid with a multitude of desirable attributes. 

One such attribute is the bound on minimum and maximum spacing. This is 

accomplished with terms in the position of A in Eq. 96. The actual construc- 

tion, however, is more readily accomplished in the finite difference form 
which is given by 


ds . 
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This represents the extended form of the earlier finite difference equation, 
Eq. 53. In the case of maximum spacing, the term 
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is considered and, within A, yield3 the contribution 
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The positive constants c and a are respectively for scaling and intensity. 
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Here, we assume the monotonically Increasing order for as displayed in Eq. 

77. If the spacing (3 i + 1 - s, ) exceeds the specified (As) , then the cor- 

responding ratio appearing in Eq. 99 exceeds unity and is further amplified by 

a which is usually about H. As a consequence, the first term in Eq. 99 gives 

the velocity of s^ a positive contribution. This is then a force to move 

in the positive direction towards s i+1 . Such a movement, of course, will shrink 

the interval (s,., - s.- ) which had exceeded the specified value (As) 

i + i i max 

Likewise, had (s ; - s,_, ) exceeded (As) , then the second terra in Eq. 99 
1 i i max 

would have caused force in the negative direction to push towards to 

thereby shorten that interval. If the intervals on both sides of exceed 

(As) , then the forces may balance each other. However, in a point itera- 
max 

tive sense, a translation in i will eventually produce an unbalanced situation 
that will successively be propagated to produce the desired result. 

In a similar manner, the minimum possible spacing is enforced with a term 
of the form 
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where c and a are the 3ame sort of positive constants and (As) . is the spe- 

min 

cified lower bound on the spacing. As above, we get the contribution 
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to the velocity of s^. Because of the minus sign, the terms are now inter- 
changed. When (3 1+1 - s^ falls below (As)^, the velocity contribution in 
the second term becomes important and pushes 3^ in the negative direction 
towards s^ and thus enlarges (a i+1 - s^. When (s^ - 3^_ ^ ) falls below 
(A 3 ) m in, a significant contribution comes from the first term to push Sj in 
the positive direction towards in order to enlarge (s^ - s^_.j). These 

terms are due to Winkler, Mihalas, and Norman [211] and were called floors and 
ceiling 3 on the spacing. 


- 38 - 



Eccentricity 


A certain amount of smoothness in the gradations of the grid can be con- 
trolled by a force which tries to equally distribute the ratios of eccentric- 
ity in adjacent spacings over the space of indices £ - i. This equidistribu- 
tion over £ rather than s is represented by the function A in Eqs. 96-97. At 
a point j, the spacing eccentricity is given by 

c ■ - h 

J Sj - s jM 0 02) 

and is not entirely suitable for immediate use because of the directional bias 
represented by a separation of j-1 and J + 1 values between the denominator and 
the numerator. Otherwise, it would have been a good control to provide smooth 
gradations in spacing. Fortunately, the eccentricity ratio c ^/e. removes 

the directional bias and provides the same control. That ratio term is given 

by 


IV 

‘j + 1/2 


(s j*t - 3 i>* 


■ 3 j*2 - S J.1 ,l! j - 3 j-,' 


(103) 


where c is a positive scaling constant. Upon equidistribution over the i, we 
would get the geometric mean 

E J = /e J-1 E J + 1 (1 0*1) 

to provide the growth rate. In such a situation, the distribution would be 
determined solely from the end conditions which would also require an extra 
point at each end. The extra points are needed for end-point spacings. 

When the eccentricity ratio term of Eq. 103 i 3 inserted into the parabol- 
ic form of Eq. 97, the force on the velocity of 3l comes from the difference 
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~ 3 i^ _ K - 3 i-i )2 

^ " 9 i + 1 ^ 3 i ~ 3 i-i J ^ 3 i + i " ®i^ 3 i-i ‘ s i-2 J 


Now suppose that In a point iterative process, s i were to miss this desired 
equidistri button by getting too close to s^.,. Then the Interval length (Sj - 
s i -1 ) ^ Eq. 105 would be too small. The effect would be to enlarge the first 
term and shrink the second. That enlargement can be substantial since the 
interval smallness then appears in a denominator. In fact, this is a control 
by means of a singularity; namely, a pole. As a consequence, the first term 
contributes a positive contribution to force s^ towards s^ + ^ thus enlarging 
the offending small interval (s^ - s^). Likewise, the second term provides 
a restoring force when the point s i is too close to s i+1 . This is just the 
opposite situation. The form of smoothness here was developed by Winkler, 
Mihalas, and Norman [24]. 


METRIC NOTATION 


In the discussion of the various adaptive methods in higher dimensions, 
the fundamental properties of each object can be more clearly and concisely 
stated with the use of metric notation than without it. As a consequence, we 
shall describe just enougn of it for our purposes. In two dimensions, we as- 
sociate 1 with E and 2 with n- The differential element of arc length ds is 
then given in a plane by 


ds 2 = g n dS 2 + 2g 12 d£dn + g 22 d n 2 


( 106 ) 


where 
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(107) 


T 

and r = (x,y) is the position vector. This form of ds 2 can be readily veri- 
fied by a chain rule expansion of dx 2 + dy 2 . The matrix of elements g A j i 3 
called the metric because it defines the distance measurements with respect to 
the coordinates (£,n). In a straightforward manner, the determinant 
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g - det(g 1 j ) 


( 108 ) 


is also seen to be the square of the Jacobian 

J - x_y - x y (109) 

For grids, this gives the cell sizes in the physical (x,y)-space that corres- 
pond with the fixed uniformly sized cells in the computational or logical 
space (E,n). When the Jacobian is nonvanishing, the inverse of the metric 
matrix (g^) can be computed and has elements that are typically written in 
the superscripted format g 1J . In continuation, the metric notation is equally 
valid in higher dimensions. For surfaces in three dimensions, we merely ex- 
tend Eq. 107 by including terras for z that come from the dot products between 
combinations of the natural curvewise tangents r and where r « (x,y,z). 

For volumes in three dimensions, we just associate 3 with a third variable £ 
and then repeat the above discussion where, of course, the expressions are 
slightly longer but the meaning remains the same. For a more compact nota- 
tion, we shall sometimes use ^ ^ f>or h, C and X 1 , X 2 , X^ for x, y, 

z. Further discussion can be found in virtually any text on differential 
geometry. A development specifically aimed at grid generation can be found in 
either of the monographs by Eiseman [25] and Warsi [26]. 

CURVE BY CURVE METHODS 

With a basic understanding established, we next proceed to consider high- 
er dimensions. The most direct extension into higher dimensions is to apply 
adaptivity on a curve by curve basis and to cycle through one or more coordi- 
nate directions. The methods of this description are called "alternating 
direction adaptive methods." From a geometric viewpoint, these methods oper- 
ate by specifying the diagonal part of the metric tensor since each such di- 
agonal entry is inversely proportional to a specified weight. This funda- 
mental geometric framework provides a unified setting for all such methods and 
was presented by Eiseman [27]. The unity here derives from the fact that no 
matter how the metric specification is accomplished, the result must be the 
same up to the distinctive errors of approximation. The completeness in the 
specif ication comes from the curvature of space. This is most readily seen in 
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two dimensions where Gaussian curvature provides a single relationship between 
all three metric coefficients g 1 1 , g 22 and 812 * That raean 3 that the coordi- 
nates are analytically determined when and g 22 are specified as positive 
functions that may even involve 8 i 2 * In other words, the two available de- 
grees of freedom are employed. Likewise, in three dimensions, the three 
available degrees of freedom are consumed. The typical specifications include 
gradient magnitudes, various forms for curvature (normal, geodesic, mean, 
second derivatives), cell properties (eccentricity, aspect ratio, lengths, 
area or volume), and the attraction to prescribed distributions (uniform, 
arbitrary, orthogonal alignment, or previous locations). Such specifications 
usually appear in a weight function. The most commonly utilized form is the 
linear one given by Eq. 1. Altogether, there is only one weight function em- 
ployed for each degree of freedom consumed. 

Metric Specification 


The main distinguishing feature from the previously examined isolated 
curves is that the curves considered here appear as part of a higher dimen- 
sional coordinate system. To appropriately account for this fact, a direc- 
tional index must first be attached to the earlier equidistribution statements 
for curves. In particular, the basic differential statement of Eq. 2 becomes 

w . ds . » c . d£ , (110) 

ii l i 

where i is now the coordinate directional index. For each given curve, c^ is 
a constant as earlier. However, as we go from curve to curve within the fam- 
ily for a given direction, that constant becomes a function of the transverse 
variables which govern this progression. As a consequence, the condition on 
Cj is that it vanish under differentiation with respect to the curve variable 

V 

Along the curve, the incremental distance measurement ds^ is just a spe- 
cial case of the measurement along arbitrary curves as indicated in Eq. 106 
and its discussed extensions. Quite simply, with the exception of dE^, all of 
the remaining differentials d£ k for k - i vanish. This leads to the form 
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(Ill) 


da. • dSj 


Upon substitution into the equidistribution statement of Eq. 110, a cancella- 
tion of d£ leads to the basic metric statement 
1 


11 



(112) 


which is the direct assertion that a prescribed weight w^ is equivalent to a 
prescribed metric coefficient g^. 

When the curve constant c^ is isolated from the metric and the weight 
w py solving Eq. 112 for c 2 , a differentiation in then removes that con- 
stant and yields the partial differential equation 


9 *ii + 2 3w i 
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(113) 


This corresponds to the earlier ordinary differential equation of Eq. 23. As 
in the case of Eq. 23, the data from the constant has been effectively trans- 
ferred into boundary conditions. Upon substitution for g ^ in Eq. 113# the 
partial differential equation is readily observed to be second-order. To be 
explicit on this matter, we consider the case of the first equation in two- 
dimensional Euclidean space. There, the metric from Eq. 107 is just g^ * 

x 2 + y 2 and results in 

5 E 
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The equation for the q-direction ha3 the same form and is obtained similarly. 
Other forms cover curves on surfaces and higher spatial dimension. All of 
these come from Eq. 113. In comparison with the earlier case of an isolated 
curve, the equations here are for coordinate curves, and as a consequence, 
reflect their embedding within the coordinating system. 


The Historical Development of Pure Curve by Curve Methods 

To examine the historical roots of alternating direction adaptivity, we 
return to one dimension and then expand to higher dimensions and to multiple 
directions. Some of the earlier studies in one dimension were developed by 
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Winkler [ 28 ], Ablow and Schecter [29], and White [30,31]. Winkler [28] con- 
sidered grid point movement directly in physical 3pace in response to gra- 
dients in the monitor surface. All attributes for clustering and most for 
grid structure were given in the form of a linear weight. By contrast, White 
[ 30 , 31 ] developed his grid directly on the solution curve and thereby effec- 
tively used a weight of unity. When the arc length was expressed as an ex- 
plicit function of physical space, he obtained the appropriate weight function 
for arc length. In extrapolating from here, he called the weights monitor 
functions. Although the term may be appropriately descriptive when everything 
is combined under a single integral, it is otherwise deficient. A more basic 
consolidation of the data and one that applies regardless of dimensionality is 
the concept of monitor surface as discussed in the introduction. Ablow and 
Schecter [ 29 ] preceded White and considered a linear weight with curvature 
that was applied relative to the arc length of the solution curve. 

In the next stage of development, Dwyer et al . [17,32~35] considered the 
process of adapting the points along each coordinate curve in a fixed direc- 
tion. In contrast to Winkler, White and Ablow and Schecter, he considered 
weight functions that depended upon positions in physical space. This was 
executed in a noniterative fashion by employing the backward global integral 
statement of Eq. 6 along each coordinate curve in the physical region. Again, 
linear weights were employed. In terms of Eq. 1 , he considered up to two 
masses consisting of magnitudes for first and second derivatives of the moni- 
tor surface. These, however, provided only approximate gradient and curvature 
data: the variations along transverse coordinate curves were ignored. An 

advantage that evolved from the choice of linear spatially dependent weights 
was the capability to more rationally define the coefficients in the linear 
weights. With the transformation established by the global integrals, Dwyer 
[ 17 ] noted that the fractional contribution of each mass was merely a ratio of 
that mass integral to the total mass integral. Each integral, of course, was 
taken over the entire curve. Since each mass integral contains the associated 
coefficients, he was able to specify the fractions and solve for the coeffi- 
cients. This was done for Up to two masses. The general case was discussed 
earlier herein. As a consequence, each specified fraction resulted in a 
weight coefficient that dynamically adjusted from curve to curve. 

In progressing from a single direction to multiple directions both Ablow 
[ 36 ] and Gnoffo [37] perfomred bidirectional studies. Ablow considered the 
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arc 


length along the curves on a monitor surface and proceeded to solve the 
equations 


Vee * Vee * Vee ■ 0 


Vee * Vee * Vee 


which come from Eq. 113 with constant weights and with z = z(x,y). He em- 
ployed an ADI procedure for the solution. As a consequence, he automatically 
clustered to high gradients. In contrast, Gnoffo viewed the monitor surface 
only from physical space, and thus had to explicitly use gradient information. 
Like Dwyer [33], he neglected transverse variations and considered derivative 
magnitudes along the given coordinate curve in physical space. He did not, 
however, consider derivatives beyond first-order. In the execution phase, the 
weight was viewed as a function of the grid point index (equivalently £ or n) 
and the forward global Integral statement of Eq. 15 was employed. That re- 
quired iteration as discussed earlier. In the trapezoidal quadrature rule for 
the integral, each increment in £ and n was unity and the result was a simple 
sum of reciprocal weights as given by Eq. 18. This is in contrast to the non- 
iterative statement where the quadrature is a sum of products between midpoint 
weight values and the corresponding interval arc lengths from Eq. 8. As a 
matter of terminology, he called this approach a spring analogy. Because of 
the iterative nature, however, the spring constants are not really constants 
3 inee they must also change. At best, they may be viewed then as nonlinear 
springs. 

In a more general study, Ei3eman [27] consolidated and extended the pre- 
vious work and presented a mathematical foundation for all such curve by curve 
methods. To be descriptive, these methods were then called alternating direc- 
tion adaptive methods. The mathematical foundation was mentioned earlier. In 
brief terms, the known curvature of a region implies that metric specifica- 
tions along each coordiante direction (Eq. 112) are enough to completely de- 
termine the metric which In turn can be employed to generate the grid by line 
integration. In the context of orthogonal grid generation, details on the 
line integration for coordinate positions can be found in Wars! and Thompson 
[38] and in Eiseman [21]. 

With the previous work separated between surface grids without weights 
and planar grids with weights, the more unified approach is to simultaneously 
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consider surfaces and weights. In the consolidated form, Eiseman [27] stated 
a preference for generating grids on monitor surfaces while using weights for 
the resolution of surface properties; albeit, the form applies equally well 
to viewing the surface from the physical region below it. The stated prefer- 
ence quite naturally derives from the fact that the accurate representation of 
a surface is more readily apparent in the surface grid than in the correspond- 
ing projected grid in the physical region. The reason is that the viewpoint 
is the location of arc length measurement which is more naturally taken and 
accurately controlled on the surface. In simple terms, it is reasonable to 
generate the grid directly upon the very object that must be given a good 
representation. 

With the objective of providing a good representation for the monitor 
surface, the geometric parameters for the surface must be used along with the 
parameters which govern the grid quality in the sense of good structure. The 
basic surface properties are the bends or folds in the surface and are the 
surface boundaries which may also be bent in some manner. The basic measure- 
ment for surface bending in a given direction is the normal curvature. The 
measurement for the boundaries is the geodesic curvature. Each curvature de- 
tects only the desired property. Using the normal curvature in the direction 
of the current coordinate curve, Eiseman [27] observed the desired clustering 
effects for surface folds. In addition, the use of normal curvature was seen 
to provide some alignment of coordinate curves with folds in the surface. 
Moreover, the formulation used the general linear weight of Eq. 1 where normal 
and geodesic curvatures are balanced with the unity for uniformity and a mas3 
for orthogonality attraction. The coefficient for geodesic curvature was pre- 
sented in a form that decayed upon leaving boundaries so that bent interior 
curves appearing from the iteration would not unnaturally cause clustering. 

In the application, the backward global integral statement of Eq. 6 was em- 
ployed along each curve. 

The use of fractional specifications initiated by Dwyer [17] were also 
extended. In a basic sense, it was noted that the resolution of a property 
which appears over a fixed length could be depreciated by merely increasing 
the total length (or other masses). Thus a methanism was required to appro- 
priately adjust the clustering Intensity to treat the local property in the 
same manner, regardless of length. For normal curvatures we must somehow de 
tect. the likely presence of local clustering regions and the possible exces- 
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3 ive consumption of surface arc length. This was done by first forming the 
ratio R of the actual arc length of a coordinate curve on the surface to the 
minimum possible length. The latter is computed in the form of a Euclidean 
distance using the arc length of the projected curve and the change in alti- 
tude between endpoints. Upon application, a factor of the form tanh[D(R-1)] 
was applied to a constant fraction determined in the original manner. The 
fraction is now seen as a maximum possible fraction that would be employed in 
the most severe case within the family of all coordinate curves in a given 
direction. The constant D then gives the rapidity for which the specified 
maximum fraction is approached. More details are provided in Eiseman [27]. 

The extension of fractional specifications to include any number of masses in 
Eq. 1 and to consider what happens when distinct masses appear on the same 
interval are all discussed in the review by Eiseman [8] and in some detail 
here in the sections about Eqs. 58-72. 

In several subsequent studies, Nakahashi and Deiwert [20,39,^0] added a 
few more items of interest to the development of alternating direction adap 
tivity, and in addition, presented some rather good examples of adaptive 
simulation in aeronautics. The main contributed item is their incorporation 
of an orthogonality control outside of the weights. Given the arc length 
locations ^ of Eq. 78 corresponding to an orthogonal alignment with the 
previous curve, they added a term of the form D^(s^ — r^) to the right hand 
side of the tridiagonal system for the direct finite difference statement of 
Eq. 29. They arrived at Eq. 83 where an increase in causes the diagonal 
dominance to grow and thus forces s i to approach r^. A variant is to balance 
the orthogonal locations with the straight line continuations through the 
previous curve. The purpose wa 3 to choose r^ so that the continued transverse 
curves would be smoother. In the applications, actually varied inversely 
with respect to the distance from the previous curve at each i. This streng- 
thened the attraction for closely spaced curves. 

In three dimensions, there is another similar arc length location t i and 
another term E^(s^ - t^). There is added to the diagonal and D^r^ + 

g t . forms the right-hand side. In continuing the analogy with springs, they 
attributed the orthogonality control to torsion springs rather than to diag- 
onal dominance. Although this control might at first appear to make the 
method distinctive, the fundamental fact still remains that a metric rela- 
tionship is being specified along curves. In fact, as we saw in Eqs. 84-87, 
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the diagonal dominance control can be interpreted as a mass for orthogonality 
attraction as discussed in regard to the general linear weight of Eq. 1 . When 
Dj is chosen as in Eq. 87, the quadrature, Eq. 8, of the backward global 
integral statement can also be employed because unique half-point weights are 
then determined. 

Rather than a consistent use of linear weights, Nakahashi and Deiwert 
[40] considered weights where some of the specified constants appeared nonlin- 
early. This arose mainly because the minimum and maximum spacing could be 
specified by means of a single algebraic formula over the curve. This was 
distinct from their use of ratios as discussed in Eq3. 73~76. One constant 
was an exponent and had to be determined by iteration. By contrast, Winkler , 
Mihalas and Norman [2*0 gave upper and lower bounds upon the spacing within 
the context of the general linear weight of Eq. 1 and did not require an 
iterative determination of constants. This was possible because the spacing 
was only required to become close to minimum and maximum spacing along the 
curve rather than matching it precisely. The development is presented here 
about Eqs. 97-101. 

It should be noted that the objective in setting bounds upon the spacing 
as represented by Winkler, Mihalas and Norman [24,41] and by Nakahashi and 
Deiwert [40] is essentially the same objective as specifying the fractional 
amount of quantities as represented by Dwyer [17] and Ei3eman [8,27]. Both 
prescriptions merely attempt to set limits upon the finite distribution of 
points. In comparison, the spacing bounds are attractive because of their 
direct attachment the the actual spacing while the specification of the 
fractions are attractive because of their flexibility. With the fractions, 
the constraints upon spacing can be more effectively balanced against the 
other attributes which quite naturally enter into the same linear format. 
Moreover, those other attributes can also be precisely separated at the same 
time as the spacing requirements. This extensibility is not readily apparent 
in the format of spacing bounds. To consider it, the form presented by 
Winkler, Mihalas and Norman [24,41] would be preferred because it fits nicely 
into the context at the general parabolic evolution equation, Eq. 97. 

Unlike the study of Eiseman [27], Nakahashi and Deiwert [20,39,40] did 
not acknowledge any defect in the application of the method. The only in- 
dication appears indirectly when they state that some corrective action is re- 
quired when the orthogonally aligned arc length locations or tj fail to be 



monotonic. Although the action was simply executed, the real problem was 
covered up. The real problem comes from the errors incurred in the numerous 
piecewise linear approximations to curves and their parameter i zat ions . Vari- 
ous forms of the problem were illustrated by Eiseraan [27] along with appro- 
priate explicit corrective actions. In cases with rapid but not excessive 
variations, such actions, however, are usually not required. 

Curve by Curve Adaptation with Derivative Smoothness 

While the straight curve by curve adaptation yields good results in many 
circumstances, there is a significant underlying limitation on the weights. 
Namely, the weights cannot be too severe or else the procedure will collapse. 
This has been observed and while corrective action can be inserted directly 
into the process, such action is rather detailed and technical. 

In a further study, Eiseman [M2] found that a better course of action is 
to redefine the directional sweeps by splitting them into two phases: the ac- 

tive phase and the passive phase. In the active phase we just have the origi- 
nal curve by curve strategy in the current direction. This contains the fun- 
damental adaptive forces. In the passive phase, a low pass filter is applied 
to remove any wiggles or abrupt changes in spacing caused by the active phase 
but to leave intact the basic results of the intended action. This produces a 
smooth grid in the sense of derivative continuity. A3 a consequence, continu- 
ous numerical derivatives are available for numerical solution algorithms and 
for the application of controls in successive sweeps. Such controls include 
the use of orthogonality and curvature in the weights. As a practical matter, 
it has been observed that the splitting of sweeps into this pr edi ct or- correc- 
tor format of active and passive phases has resulted in considerably enhanced 
stability and a much larger range of severity in the choice of weights. 

The simplest form for the passive phase is given by the direct action of 
a Laplace filter upon the grid point locations. While such action by itself 
may not be appropriate for the generation of a nonsingular grid, it is cer- 
tainly suitable for the stated purpose of establishing derivative continuity. 
The distinction is that the filter is applied at mo3t a few times in each di- 
rectional sweep rather than being driven towards convergence as is the case 
when the solution to a system of grid generation equations i3 sought. The 
Laplace filter employed by Eiseman [M2] was given by the simple Gauss-Seidel 
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At the boundaries, this formula was restricted to provide filtering in only 
tangential directions. In the application, a three-dimensional monitor sur- 
face was defined in four-dimensional Euclidean space by means of the vector 
(x,y,z,u) where u is considered to be a function of (x,y,z). To obtain the 
physical space projection, we merely replace u with a constant which is usual- 
ly 0. In a test of the basic movement, u was taken to have a severe variation 
across two intersecting ellipsoids that also intersected the boundaries of a 
Cartesian box. From an initial surface grid with a Cartesian grid projection, 
an equal arc length grid was rapidly generated on the surface, and according- 
ly, a smooth gradient clustered grid resulted in the (x,y ,z )-projection 
(x,y ,z,0) . 


The Conformal Measure of Smoothness: Basic Elliptic Grid Generation 


While the passive phase of each sweep provided smoothness in the sense of 
derivative continuity, another measure of smoothness is provided by an at- 
traction to conformal conditions. This latter measure is more demanding and 
is thus clearly stronger. 

The conformal measure of smoothness is most readily observed in a varia- 
tional setting [8]. In two dimensions, we simply minimize the amount by which 
the Cauohy-Riemann conditions fail to be satisfied when generally incompatible 
boundary conditions are employed. At a point In the field, each of the 
Cauchy- R iemann equations has a residual that in the general circumstance 
deviates from zero. A measure of nonsati3faction is clearly given by the sum 
of squared residuals. For a mapping from (x,y) to (£,n)i the smallest loss 
over the whole field is then obtained when the integral 

T f (U - n ) 3 ♦ U + n )* ]dxdy (117) 

c ; x y y x 

i 3 minimized. Although the mapping direction in the formulation is opposite 
to that of an eventual application, there is a good reason: namely, that the 
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control is relative to whole coordinate curves that are associated with con- 
stant values of £ and n rather than to horizontal or vertical lines for con- 
stant values of x and y. The coordinate curves in the physical space (x,y) 
are generally not straight lines. More discussion of this choice can be found 
in a variety of sources [6-9,12]. 

When the integrand in Eq. 117 is expanded and regrouped, we get 

I = J[U 2 + f; 2 ) + (n 2 + n 2 ) + 2(5n - S n )]dxdy ( 11 8 ) 

c J x y x y k y x x y 

where we can identify the newly grouped terms. From Eqs. 107 and 109, we get 

I c = / [g 11 * g 22 ~ -jjdxdy x (119) 

Since dxdy = Jd^dn, the integral simplifies to be 

I r = j (g 11 + g 22 )dxdy - 2 j d£dn (120) 


where now the last term is a constant because the domain of curvilinear var- 
iables U,n) has a constant area: typically, this is just a rectangle or some 

collection of them. As a consequence, conformal smoothness is obtained when 
the constant part of Eq. 120 is dropped and the remaining integral 

1^ " j (g 11 + g 22 )dxdy (1 21 ) 


is minimized. In this form, the extension into three dimensions is clear: we 

just add a term g^3 to the integrand and integrate over volumes. As in Eq. 

50, the Euler equations for the minimization of I g are obtained by applying 
the general operator 


3x 


i 


»(x|> 


( 122 ) 


to the integrand of Eq. 121 and setting the result equal to zero. For each 
equation i, the repeated index j in the operator is assumed to be summed from 
1 to the number of spatial dimensions. This is just a commonly employed sum- 
mation convention. For Eq. 121, we get the Laplace system 
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( 123 ) 


that was considered by Winslow [43]. Alternatively, the same system could 
have been derived by transforming the integral I 3 to curvilinear variables 
5^ * 5 and ^2 " n and then em Pl°ying the general operator 


__3_ _ _3_ 3 

3 5 • 3x . 3(5 ) 
i J i x 


(124) 


to the adjusted integrand. This is just the 
fact that the Laplace formulation of Eq. 123 
variables; Thompson, Thames and Mastin [44] 
zat ion 


parallel to Eq. 51 . Using the 
is directly for the curvilinear 
considered the Poisson generali- 


V 2 5 = p 

V ? n - Q 


(125) 


where the forcing terms P and Q provided separate and active controls on the 
respective families of coordinate curves. Because of the general effective- 
ness of this grid generation system, it became thoroughly developed and widely 
used. The effectiveness arose primarily from the treatment of whole coordi- 
nate curves relative to the conformal attraction inherent in Eq. 123. To ap- 
ply the Poisson system, Eq. 125, the dependent and independent variables must 
be Interchanged. The result is given by 


11 ^ 12 22 „ „ 
g x + 2g x + g x + x P + x Q * 0 

55 5n nn 5 n 

11 12 22 
g y rr + 2g y - + g y„„ + y r P + y Q - o 
55 5n nn 5 n 


( 126 ) 


In the general case, we have 

g lj r + P.r HH (12 7) 

hh 1 h 

where r is the position vector in physical space, is the i~direction forc- 
ing function and H is the mean curvature should the application be for a 
surface with unit normal N. In a Euclidean 3pace, H vanishes and the right- 
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hand side then also vanishes. 

Upon examining Eq. 126, Middlecoff and Thomas [45] noticed that the forc- 
ing functions would have to be appropriately scaled to account for potentially 
disparate variations in the inverse metric elements g*t As a consequence, 
they redefined the forcing terms to include a factor of g 11 for each direction 
i. In terms of the general form of Eq. 127, we get 


g lj (r r r + 5 . .7 . r ) » HN (128) 

5 1J 1 L 

1 J 1 

where the adjusted forcing functions are given by 



and the <5 in the sum on i and j is the Kronecker symbol that is unity when i 
ij 

= j and vanishes otherwise. For more details on the general development of 
this Poisson system, the reader is referred to the review articles by Thomp- 
son, Warsi and Mastin [6], Thompson [7], Eiseman [8], and Eiseman and Erle- 
bacher [9] as well as the text by Thompson, Warsi and Mastin [46]. 


Curve by Curve Adaptation with Conformal Smoothness 


Recognizing the basic need for smoothness, Anderson and Steinbrenner 
[47,48] brought the process of equidistributing weights along coordinate 
curves into the format of the Poisson system of Eq. 125. Their development 
was motivated by the previous work of Middlecoff and Thomas [45] who employed 
the formulation of Eq. 128 in the planar form 




♦ r E ) ■ 26 12% * e n ( 


nn 


Vr ) 
n 


( 130 ) 


with 
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for the purpose of converting boundary distributions into forcing functions 
thereon. While Middlecoff and Thomas interpolated the forcing functions into 
the field so that boundary distributions could be propagated into the interior 
of the field, Anderson and Steinbrenner viewed each curve as if it were such a 
boundary curve. In the boundary application, the assumption of orthogonality 
and vanishing transverse curvature was employed to get the differential equa- 
tion for an equidistributed weight. In particular, with g 12 - r *r^ * 0 and 

r =0, the eauidistr i bution statement of Eqs. 2H-25 is obtained from Eq. 1 30 
nn 

as 



where s is the curve arc length and 



(1 32 ) 


033 ) 


is the relationship bewteen the forcing function $ and the weight function w. 
In the more general interior application, the same sort of equidistribution 
statement was established without the previous assumptions. The main dis- 
tinction is the addition to $ of a term for orthogonality and terms for curva- 
ture. These orthogonality and curvature terms represent the attraction to the 
conformal measure of smoothness. Without forces, they yield the desired 
smoothness. With forces a deviation is obtained. 

In the general two-dimensional case represented by Eq. 130, our task is 
to obtain a reduction into the form of Eqs. 132-133 for each direction. To 
start, we shall rewrite Eq. 130 in a form where the basic constituent vectors 
are the unit tangent vectors along coordinate curves and the unit normal vec- 
tors perpendicular to coordinate curves. From the metric coefficient defini- 
tion in Eq. 107, the appropriate normalization factors for the natural tangent 

directions r and r are obtained directly. With t. and denoting the res- 
E n 1 d 

pective unit tangents, we have 
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( 13 *») 


where s - /g^ and u^ - /g^ are the respective arc length derivatives. For 
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Eq. 130, second derivatives are needed and can be obtained by derivatives of 
Eq. 134. To get them, however, the unit tangents must be differentiated. 
This is accomplished with the first of the Frenet formulas (cf. [49]) which 
are given by 
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for the two respective coordinate directions. The unit normal vectors and 
curve curvatures are given by and and by k 1 and k 2 . By using the chain 
rule, the arc length dif f erentiation in Eq. 135 is readily converted into dif- 
ferention with respect to Z and n. With the result inserted into the deriva- 
tives of Eq. 134, we arrive at 
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Upon substitution of the first derivatives in Eq. 134 and the second deriva- 
tives in Eq. 136, the grid generation equations of Eq. 130 become 
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where the terms have been grouped according to the tangent and normal direc- 
tions. An immediate observation from the grouping is that the equidistribu- 
tion process along coordinate curves is now separated by direction: namely, 

that the repr esentati ve form appears only in the coefficients of and x . 
To isolate the forcing functions, only a dot product with the respective 
normal vectors and is required. To remove ¥, a dot product with 
yields 
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As a matter of interpretation, the deviation from an equidistribution of <j> is 
represented by the second term of Eq. 1 ^0. Xt contains an orthogonality part 
due to and *n^ and curvature parts due to and k 2 for the coordinate 

curves in the two distinct directions. When the coordinates are orthogonal, 
that deviation reduces to On balance, the force towards an equi- 

distribution of <j> must overcome the forces of conformality that are repre- 
sented by orthogonality and curvature. A similar development and balance oc- 
curs for 

In the adaptive context, when the local force is sufficiently strong, the 
equidistribution force overpowers the smoothness conditions to become dominant 
and, thereby, to provide the desired equidistribution of the weight along each 
curve. The equidistribution is more localized than the previous derivative 
continuous measure of smoothness. This occurs because the equidistribution is 
essentially cut off unless it is sufficiently strong. The value of 3uch a 
cutoff is that intense local clusters can be formed primarily from curves that 
are not too far from the given locality. In contrast, the actual equidistri- 
bution of weights along curves adjusts all points along curves, and thus, 
tends to globally propagate adjustments to intense local requirements. This 
tendency can of course be limited with the explicit use of orthogonality and 
curvature attraction. 

The Variational Form of Curve by Curve Adaptation 

In the prior section, the curve by curve controls were Injected into the 
two-dimensional Poisson system of Eq. 125 in a somewhat ad hoc manner. This 
was done by first casting the system into a form, Eq. 137, that explicitly 
displayed the appropriate elements of equidistribution and then by making the 
ad hoc judgment that those elements were dominant. The judgment occurred at 
the stage represented by Eq. 1 40. While the same development could be con- 
tinued into three dimensions, the added detail in getting to the stage of 
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judgment, not to mention the judgment itself, la much more complex. 

In order to simplify the judgment and to gain a higher degree of general- 
ity, we shall develop a variational formulation based upon the fundamental 
statements of equi distr i bution and conformal smoothness. Upon observing the 
Integrand of the smoothness integral I of Eq. 121, it is noted that the i-th 
coordinate direction is represented by the inverse metric term g^. In a geo- 
metric sense, this represents the spacing between curves as we march along in 
the i-direction. To control that spacing, we consider a weight function w^ 
which would be completely effective if g 11 /w^ were constant. For orthogonal 
systems, the basic metric equidistri bution statement of Eq. 112 is obtained 

for AT . On the overall basis where the various terms in the integrand must 

i 

compete with each other, we are led to the integral 
„ 1 1 22 

I = | (£ — + ® — )dxdy 0*41) 


The associated Euler equations derived from the operator of Eq. 122 are given 
by 


( 1 * 42 ) 
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With the exception of the g 1 ^ terms, this matches the equations of Anderson 

i p 

and Steinbrenner [47,48] that were discussed in the prior section. The g 
terms, here, represent the transverse variations of the weights for each 
direction. As in the prior case, the inclusion within the Poisson format is 
clearly quite simple. When w. = w 2 » w, the conformal forces are scaled by 
1/w in the variational integral of Eq. 1 41 . Given the conformal tendency to 
maintain orthogonality, the uniform scaling of /w along each coordinate In- 
tuitively suggests an equidistribution of w with respect to volume elements. 
This equidistribution was analytically verified by Anderson [50] who was moti- 
vated by the diffusive form proposed by Winslow [51] rather than the above in- 
tuition. The diffusive form comes from the use of w * 1/D in Eqs. 141-142. 
With the prescribed diffusion function D, the equations proposed by Anderson 
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[50] are obtained by setting w 1 - w 2 « 1 /D in Eq. 132. Here, the logarithmic 
derivative trades a reciprocal for a minus sign so that the form of the forc- 
ing terms in Eq. 1 42 now appear with minus signs and a replacement of w 1 and 
Wp by D. 

Unlike the development of Anderson and Steinbrenner , the extension into 
three dimensions is easily formulated and justified. From the integral 
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we get the Euler equations 
7 2 r = g 


.. . ..i iVs. ,.j ]Vs 

w 8 W, 6 W, 


1 


1 


1 


V 2 n - g 


2, (U 2>5 22 ( “2>n „ _23 


2 


+ g 

w„ ° w. 


(133) 


(133) 


(w.) r „ (w ) (w ) 

?’t • s 3 ’ -5^ * 8 32 T 2 - 3 * s” ~^ JL 
3 3 3 

In a more compact form, the two- and three-dimensional equations can be stated 
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Upon returning to the basic statement of Eq. 112 and noting that the de- 
rived equt distribution was for each /w7 rather than Wj , we are motivated to 
replace each w, by its square and repeat the argument. Even more generally, 
it is just as easy to consider a replacement by w” for any positive constant 
Cl. When the replacement is made, we can simply carry w^ through all of the 
above steps in an undisturbed manner and then perform the final differentia- 
tion at the end. Altogether, we obtain the general forcing terms 


P =ag U 0 3 

i w . 

l 

that can be employed in Eq. 127 which also i3 equally valid on surfaces. As 
originally mentioned, when a = 2, the equidistribution forces conform most di 
rectly to the basic metric equidistribution statement of Eq. 112. This might 
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actually be the most intuitively plausible application when the linear weight 
of Eq. 1 is employed. Then, at least, there is a more direct correlation to 
the previous curve by curve studies that use linear weights. 

With the forcing terms of Eq. 1H6 in the elliptic operator from Eq. 127, 
a concise evolutionary parabolic partial differential equation can be deduced, 
In a basic sense, it is the natural extension of the one-dimensional form 
presented earlier in Eq. 52 and is given by 


a . . 'Vs. 

■ 8 lJ|r s.s. + “ ^ V 1 ' HN 

1 J 1 1 


(1 U7) 


whe~e linear weights retain the original meaning when a - 2. The further ex- 
tension in the spirit of Eq. 96 is also possible. For instance, the attrac- 
tion to a fixed prescribed grid q can be inserted a3 a term within the format 
of 3 in Eq. 96. As in Eq. 9 1 *, the grid r is attracted to the fixed q with the 
term D(r - q). The evolutionary equation then becomes 

. .. 'Vs. 

K ~ - gW (r + a J r_ ) - HN - D(r - <]) 0«S) 
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FINITE VOLUME METHODS 


Finite volume methods are methods where the grid point motion is based 
upon the volume elements between grid points. The basic format was first es- 
tablished in one dimension with the mean-value relaxation statement of equi- 
distribution. This led to the center of mas form of Eqs. 32~33 and the move- 
ment form in Eqs. 3*1-35. The second format was correspondingly established in 
the variational setting with Eq. 36. 

Direct Mean-Value Relaxation 


For simplicity, only two dimensions will be considered since the funda- 
mental pattern is established therein. A two-dimensional stencil centered 
about a grid point rv , is displayed in Figure 1. The local motion of r^j comes 
from the weighted volume elements which are defined by the indicated triangles 
containing as a vertex. In each quadrant k, a barycenter b k and a weight 
w k is obtained as the average of positions and values respectively over the 
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triangle vertices. With the weight w k considered to be uniformly distributed 
over the corresponding triangle of area A k , the contribution at the barycenter 

b, is just w,A. . From the weighted triangle areas, the direct extension of Eq. 

K K K 

32 is the center of mass for the four quadrants and is given by 


U 



Figure 1 : Finite volume stencil in two dimensions 
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to determine the new position for r 


The actual movement is given by 
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and is employed in a point iterative cycle. This represents the most direct 
statement of mean-value relaxation. In the movement form of Eq. 150, the ef- 
fective origin for the computation has been shifted to the stencil center r^j. 

Like the earlier Laplace system of Eq. 123, the motion represents an at- 
traction to conformal conditions when each w k is unity. In terms of the lin- 
ear weight form of Eq. 1 , the first term which is the number 1 is then the 
representative of the conformal attraction relative to which the other forces 
are applied. The analytical indicator for the conformal conditions comes from 
the converse of the mean-value theorem which is 3imply presented by Epstein 
[52] on pages 146-148. While the analytical argument employs the area mean 
value of functions over circular disks, the finite parallel given by Eq. 149 
with w^ = 1 is only an approximate form employed within an iterative cycle. 

The pointwise relaxation of the center of mas3 formula of Eq. 149 has 
been considered by Diaz, Kikuchi and Taylor [53], Oden, Devloo and Strouboulis 
[ 54 ], Schwartz and Connett [55], Connett, Agarwal, and Schwartz [56], Erle- 
bacher and Eiseman [57], Eiseman [58], and Erlebacher [59]. In the studies by 
Erlebacher and Eiseman, the more general application to unstructured meshe3 
was developed. 

The relaxation formula for general connectivity triangular meshes is es- 
sentially the same formula as given by Eq3. 149~150. Adjustments, however, 
must be made to accommodate any number of triangles in the determination of a 
pointwise move. This requires a careful labeling system. In our discussion, 
we shall assume a global pointwise index m for r m together with a local index 

i to label the points P = r , , . that are directly joined to r . The local 

J j a(ra,J) 

index j is attached successively to the positions as we move about p m in a 
circular manner. The range for j is taken to be from 1 to j m where J m + 1 i 3 
identified with 1 (i.e., a reduction modulo j m ). Under these index condi- 
tions, the relaxation formula becomes 
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Of the t.erni 3 represented by dot products, the first is due to the derivative 
of the m~th term in S while the second and third terms come respectively from 
derivatives of the j-th and (j+1 )-th terms that appear with distinct global 
indices in S. This occurs for each triangle j because r m appears as the cen- 
ter in the m-th term and as a surrounding term when the center is shifted to 
Pj or Here, we are assuming that r m is at least two connections away 


- 62 - 



from the boundary so that both Pj and Pj + 1 are varied and accordingly have 
corresponding sums in S. The basic coupling clearly comes from the barycen- 
tric formula of Eq. 152 and gives rise to the factors of 1/3 in Eq. 155. Us- 
ing the barycentric formula, the minimization condition of Eq. 155 reduces to 
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(156) 


which is a componentwise equality for a vector equation. The vector equation 
is obtained by just removing the dot product with in Eq. 156. The immedi- 
ate consequence is the relaxation formula stated in Eq. 151. In summary, we 
have Just shown that the center of mass formula is the general field equation 
for the variational problem stated in Eq. 153. 

Both the variational statement and the consequent field equation can be 
brought into the evolutionary context along with a mechanism for attraction to 
a given grid q. As in the curve by curve case with Eq. 147, we have an under- 
lying parabolic partial differential equation for the Euler equation. As in 
the one- di mensional evolutionary least squares statement of Eq. 57 and the as- 
sociated mean-value expression in Eq. 56, we get an extension which also in- 
corporates the attraction to the fixed q as expressed in Eq. 95. In a direct- 
ly parallel manner to the previous development, the minimization of 
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The field equation of Eq. 158, of course, will sepcialize to that of Eq. 149 
in successive stages. 

Full Mean-Value Relaxation 


In the center of mass formula of Eq. 149, the control over the grid comes 
from the specification of a single weighting quantity that is applied against 
the triangular elements of area. While this approach yields a control over 
the elemental area distribution, it does not exercise all of the available 
(jogrggs of freedom for such control. Yet one more degree of freedom is avail 
able. The full utilization of all degrees of freedom was evident in the 
transition from the Laplace system of Eq. 123 to the Poisson system of Eq. 

125. A similar transition is a reasonable expectation for the process of 
mean-value relaxation. 

The utilization of all degrees of freedom in mean-value relaxation was 
developed by Eiseman [ 60 ] . As in the Poisson system of Eq. 125 and in the al- 
ternating direction methods of adaptivity, control was established with a co- 
ordinate directional bias that provided the required separation into distinct 
degrees of freedom. In the local stencil of Figure 1, the directional bias 
was obtained by a projection of the movement vector onto the coordinate curve 
through r^j in the appropriate direction. If w is the weight for the i-diree- 
tion, then the projection is onto the curve from r ^ j to ^ j • the 

implementation, the transverse coordinate curve from j— i *”i,j + 1 used 

to divide the weights so that first and fourth quadrants would pull towards 
r. . • from r ■ while second and third quadrants would pull towards 
from r. .. On each side of the transverse coordinate curve, a center of mass 
was computed and then projected onto the appropriate segment from r^. De- 
noting the projected distances by d + and d_ for the positive and negative x 
directions from r^, the new position along the curve is given by 


d ~ 


w + d + + w_(-d_) 

w + w 

+ - 
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where 
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This i 3 the same foru as given by Eqs . 3 il_ 35 in one dimension. Since the con- 
struction is centered at r^, the sign of d gives the direction of movement: 
negative is towards r^_i j, while positive is towards r i+i > j* a 3 i m i 1 ar 

manner , a second weight w is employed for the j-direction and a similar dis- 
tance is determined along that direction. The signs of the two distances then 
determine the quadrant which contains the new position. That new position is 
determined by interpolation. 

The projected distances appearing in Eq. 159 are constructed from the 
vectors 
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which point from the current position to the nearest neighboring points 

along coordinate curves. To compute the positive projected distance d + in the 
i-direction, we must first obtain the center of gravity for the side in the 
direction of I + . This is given by 
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and reduces to 


C 


I 


+ 


r. . 

ij 




- w i A i J + 


V„j_) 


(163) 


upon employing the definition of barycenters. The projection of C 
along the positive i-direction comes from a dot product with I + /||l + || 
yields the maximum distance 
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In the original form, Eiseman [60] considered a bilinear interpolation 
that also included the diametrically opposite point to rjj. In a further 
study, Schwartz and Connnett [55] and Connett, Agarwal and Schwartz [56] ex- 
perimentally found that the method became unstable when the weights became 
sufficiently severe. While this led them to consider the earlier center of 
mass form of Eq. 149, it lead Eiseman to consider barycentric interpolation in 
place of the original bilinear interpolation. The intuitive reason was that 
the asymmetric use of diametrically opposite points would be more restrictive 
on the stability. In subsequent tests, the Intuition was confirmed: the use 

of barycentric interpolation led to stable results even when the weights had 
considerable severity. 

With the distances along the i- and j-directions determined by Eq. 159 
and denoted by d i and dj respectively, it is a simple matter to state the 
barycentric interpolation. If, for example, both dj^ and dj are positive, then 
the interpolation must be performed in the first quadrant. There, the posi- 
tively-oriented vectors I + and J + from Eq. 1 61 are employed and the interpola- 
tion is given by 
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r = r . . + d . 
ij 1J i 


J + 

nr tt + d . 
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A substitution from Eq. 1 61 provides the reduction into the barycentric form 
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where 




and 
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Since the coefficients presented by Eq. 1 67 usually lie In the unit Interval, 
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the new position from Eq. 1 66 usually lies within the triangle determined by 
the three vectors of Eq. 166 . Clearly, the same format applies to the other 
quadrants . 

Upon inspection of the full mean-value relaxation process, we note that 
the projection was accomplished with a rather strong adherence to the indi- 
vidual coordinate curves. While this is the closest to our objective of sep- 
arating the action of weights to be along the coordinate directions, a slight- 
ly weaker adherence can also be considered. The intent, here, is to get a 
simpler formulation that still retains the separation feature. The weakness 
in comparison with the more thorough version is that the center point r^j of 
the stencil in Figure 1 is removed for the projection. The consequence is 
that the algorithm may not be as robust as the stronger version; although, 
from a practical point of view, experience with the weaker formulation is seen 
to be fairly competitive. To state the weak form, an i-direction weight w is 
first employed in the mean-value movement formula of Eq. 150 and is then pro- 
jected against the i-direction secant to determine a vector increment (Ar)^ 
That is, 


(Ar ) ^ 


[(r neW 
L ij 
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where 


t . 
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r i+i j ~ r i-i ,j 


(169) 


Similarly, a j-direction weight w is used to determine (Ar)j. The new posi 
tion is then determined by the simple sum 


A r « (Ar)^ + (Ar)j 


(170) 


for the actual movement from the old r jj- 

Upon consideration of the primary motivation to separate the action of 
the weights to be along coordinate directions, we are led to a variational 
statement which is a natural extension of that in Eq. 153 Tor triangles and 
that in Eq. 36 for curves. With an impulse towards simplicity, we might first 
examine the minimization of 
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where 


n i + 1/2,J " W i>1/2J A i+1/2,j 

fl i,j+1/2 = W i , j + 1 /2 A i , j + 1 /2 

The Euler equations, however, yield the relaxation formula 
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which, unfortunately, provides a more modest separation for the coordinate 
directions. To obtain the strong separation of Eqs. 159-164, we then must 
consider the statement 


S - l (Kd + - w _ d J-j + K d + " w -0!j} (173 ' 

i ,j 

where the constructive elements about each r^j are given in Eqs. 159“l64. At 
equilibrium, we have a balance between the original movement forces and the 
boundary conditions. As a final note, here, this separation can readily be 
compared with the curve by curve approaches and extended into the evolutionary 
f or mat . 

Di rect Variational Constructions 

Rather than a direct construction of a movement formula and a subsequent 
implied variational statement, the process here will be reversed: the con- 

struction will appear in the formlation of the variational statement and the 
movement action will be the consequence. In so doing, there is the option of 
directly minimizing the defining sum rather than deriving a movement formula 
from the Euler equations. Kennon and Dulikravich [61-62] pursued the direct 
approach by using an optimization technique. The primary motivation for their 
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minimization problem came from previous variational studies [63“68] along with 
the desire to place a larger constant at the center location rjj than would 
have come from the traditional use of central differences that would have 
missed r^. This then led them to consider a finite volume stencil. 

The basic attractive forces come from distinct sums of squares that are 
balanced against each other. Of the forces, the first one pulls towards an 
equal volume distribution with terms of the form 


(a 1 -a 2 ) 2 * (a 2 -a :( ) 2 ♦ (a^-aJ 2 ♦ (a 14 -a i )- 


2 3 
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for each point r-,. The second one pulls towards an orthogonal grid with 


pO nit i i j 

terms of the form 




(175) 


where the local vectors about r^j are defined in Eq. 161. Denoting the volume 
equalization and orthogonality sums by S y and Sq respectively, Kennon and 
Dulikravich minimized the linear combination 


(1-a)S y ♦ aS Q 


(176) 


to generate or improve a grid. The parameter 0 £ o £ 1 was chosen to balance 
the two forces. The inclusion of adaptive forces comes with the inclusion of 
yet another sum. This corresponds with the basic pattern of most variational 
methods. The typical choice for an adaptive sum S a is the one over weighted 
volumes which can be assembled with terms of the form 


w(A 1 ♦ A 2 + A 3 + Ajj) 


(177) 


for each location r^. The overall sum then appears in the form 


A, S + A.S n + AS 
1 v 2 0 3 a 


( 178 ) 


which upon minimization provides a competition between the various effects. 

As in all variational constructions, various effects can be readily in- 
serted into the formulation. Here, many such effects can be developed from 


- 69 - 



the vectors of Eq. 1 61 which comprise the basic metric structure. In contrast 
with the ease of formulation, the actual use of such formulations can be 
complicated; thus, presenting additional issues. 

VARIATIONAL METHODS 

The basic pattern of developments with variational methods is to gather 
the desired attributes, form a positive pointwise measure of each such at- 
tribute, integrate the measures over the field, form a positive linear com- 
bination of the integrals, and then minimize the resultant combination. The 
minimization process typically follows the standard manipulations from the 
calculus of variations [ 1 6 3 . This leads to a system of partial differential 
equations known as Euler equations (Eqs. 50, 51, 122,12*1) that are then to be 
solved by the available numerical methods for PDE’s. The main attribute in 
most methods is the attraction to conformality and is given by the integral of 
Eq. 121 for the Laplace system of Eq. 123. This means that other attributes 
are balanced against the Laplace system for curvilinear variables that was em- 
ployed by Winslow [ 43 ], Thompson, Thames and Mastin [44] and others. The main 
adaptive attribute typically comes in the form of a weighted Jacobian so that 
the minimized integral by itself would produce an equidistribution of the 
weight over volume elements or powers of them. Along with the attributes of 
conformality and weighted Jacobians, Yanenko et al [67] included a Lagrangian 
attraction for fluid motion while Brackbill and Saltzman [ 63 ] inlcuded an or- 
thogonality attraction for an improved grid structure. The orthogonality at- 
traction came from integrals of squared cross metrics. In two dimensions, 
these were either (g^„) 2 or (8^) 2, With the motivation from the use of or- 
thogonality in Brackbill and Saltzman, Kennon and Dulikravich [61 3 pursued the 
finite volume approach discussed earlier, Saunders [68] examined a tensor pro- 
duct of B-splines, Kreis, Thames and Hassan [69] considered scaling problems, 
and Steinberg and Roache [15] developed a variant with reference grids. 

Rather than employ the attraction to the preferred Laplace system for El- 
liptic grid generation, Morice [70], Oskam and Huizing [71], and Steinberg and 
Roache [15] used attraction towards the Laplace system for locations in physi- 
cal space. While this provides conformal smoothness, this is not as robust as 
the preferred system: it does however offer some simplicity. The main in- 

surance against grid folding then falls upon the volume control or the permis- 
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sion of boundary point motion. Steinberg and Roache [15] favor volume control 
while Morice [70] and Oskam and Huizing [71] rely primarily upon boundary 
point motion. 

The idea of reference grids employed by Steinberg and Roache [15] is the 
same as that employed earlier by Steger and Chaussee [72] in their development 
of hyperbolic methods. The use, however, is more extensive in that various 
properties are extracted from the reference grids. They appear essentially in 
the form of multiplicative factors applied to the standard format of the other 
investigators. These factors are weights that give appropriate ratios between 
the current desired grid and the known reference grid. A note of caution to 
be observed in the extraction of properties is that a desired attribute might 
not be obtained when the number of available degrees of freedom is exceeded. 

A somewhat general development of the variational methods is given in the 
text by Thompson, Warsi, and Mastin [46]. This was subsequently implemented 
in a large program [731* Xn one dimension, earlier studies were given by 
Gough, Spiegel, and Toomre [74] and by Pierson and Kutler [75]. Also, a study 
of temporal smoothness was developed by Bell and Shubin [22]. In the multidx 
mensional context, a study for systematically dealing with the inherent com 
plexity of variational methods was performed by Steinberg and Roache [76,77] 
who proposed appropriate symbolic manipulation methods. 


TEMPORAL ASPECTS 

When the temporal accuracy of a simulation is to be enhanced, the evolv- 
ing grid must closely follow the trajectories of the severe disturbances in 
the solution. This is in contrast to the situation where a steady-state solu- 
tion is sought and the grid there eventually settles down into a virtually 
final configuration. The primary concern in the development of an accurate 
temporal resolution is that the severe disturbances will not escape from their 
resolution during the course of any time step in the numerical simulation. 

With this concern in mind, a number of investigators have felt that the grid 
equations should be formulated directly for grid point velocities rather than 
positions. Then, at least, the grid velocities would follow an analytical 
model for any chosen time level such as the full implicit level or the level 
halfway between explicit and fully implicit. Clearly, the use of such veloci- 
ties should be an improvement over using backward differences in time to esti 
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mate the same velocities from a grid point movement scheme that only produces 
pointwise locations from currently available data. From the viewpoint of grid 
velocities, the schemes which produce only pointwise locations can be referred 
to as static rather than dynamic (e.g., note Hyman and Naughton [78]). A typ- 
ical difference is that static methods depend upon data at only one instant of 
time while dynamic methods often depend upon data over an Interval of time. 
With the advantage of directly obtaining grid velocities, the dynamic methods 
can possess the disadvantage of a more difficult control over coordinate non- 
singularity. This is because grid point locations must ultimately be deter- 
mined, and if some velocities are too large or change too quickly, then points 
could be either overrun or artificially encircled: corresopndingly , there 

would be grid overlap or excessive skewness. 

Among the dynamic grid velocity methods, there are methods proposed by 
Winkler, Mihalas and Norman [2*1,28,41], and Bell and Shubin [22], Hindman and 
Spencer [79], Rai and Anderson [23,80,81], Greenberg [82], Piva, DiCarlo, 
Favini and Gui [ 83 ], Ghia, Ghia and Shin [84], and Harten and Hyman [85]. 
Winkler, Mihalas and Norman [ 24 , 41 ] develop a scheme based upon the equidis- 
tribution of weights over grid point indices. Nonsingularity in their one di 
mens ional context is provided primarily by creating singularity barriers which 
keep the points from interchanging positions. This comes from equidistribut- 
ing the cell eccentricity as discussed earlier with Eqs. 102-105. In addi- 
tion, they consider asymmetric time filtering to preserve resolution for the 
rapid reappearance of salient phenomena such as in shock wave reflections. 

The chief mechanism is a diffusion coefficient arising from a constructed fac- 
tor (K in Eq. 97) on the time derivative. The factor contains enough residual 
memory to slow down the diffusion of the resolution just enough to allow for 
the rapid reappearance of small structures. Otherwise, resolution would have 
to rapidly disappear and then reappear, thus causing unnecessary numerical er- 
rors because of the temporal jerkiness. In comparison. Bell and Shubin [22] 
remove temporal jerkiness by balancing the weight function equldistribution 
against the mangitudes of time derivatives in the variational format given in 
Eqs. 89 and 90. Their balancing coefficient was, however, only a constant and 
thus did not contain the residual memory as in the case of Winkler, Mihalas 
and Norman. 

In another direction, Hindman and Spencer [79] converted the Poisson sys- 
tem of Eq. 125 into a grid velocity equation by formal temporal differentia- 
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tion of the original equation. In addition, they al 30 explored the use of 
equidistributed weight functions. They found the same relationship to the 
forcing function in one dimension that Anderson and Steinbrenner [47,48] 
eventually discovered in two dimensions. 

Identifications with alternate metaphors, such as the earlier spring 
analogy, have also provided the inspiration for several methods. In this 
spirit, Rai and Anderson [23,80,81] developed an analogy to a gravitational 
potential while Greenberg [82] related the grid movement to chemical reac- 
tions. The gravitational forces decayed according ot an inverse power law for 
distances in the space of grid point indices. Force magnitudes and directions 
at each grid point came from the deviation of the error indicator from its 
average value. This was discussed earlier with Eqs. 9 1 — 93 • In a similar 
manner, the chemical reaction rates contained the adaptive data. 

An even more direct use of physically based motivation occurred in the 
somewhat parallel studies of Piva et al . [ 83 ] and Ghia, Ghia and Shin [84]. 
There, idealized two-dimensional momentum equations for viscous flow were 
transformed into diffusion equations. This was done because diffusion equa- 
tions are easier to solve. The process basically amounts to a removal of the 
convective terms which would appear when the equations are expressed in terms 
of an arbitrary time-dependent coordinate system. The two resulting equations 
are the grid movement equations which assume the Poisson format. In a sense, 
thi 3 is similar to the pursuit of Hindman and Spencer [79], although there is 

no consideration of equidistribution. 

In addition to the static methods based upon the previous solution step 
and the dynamic grid velocity methods, there are methods which impose a grid 
distribution mechanism at some implicit level without the direct determination 
of a grid velocity. This includes the methods such as that employed by White 
[ 31 ] or Mueller and Carey [ 86 ] and that employed with the moving finite-ele- 
ment method investigated by Miller and Miller [87], Miller [ 88 ], Miller [89], 
Gelinas, Doss and Miller [90], Herbst, Mitchell and Schoobie [91], Baines 
[92], and Baines and Wathen [93]. In the case of White [31] and others like 
it, the grid equation appears as a time- dependent constraint which is applied 
in an implicit coupled manner with the evolutionary physical equation. In the 
moving finite-element method, the coupling and the grid evolution comes 
directly out of the formal finite-element process when it is directly extended 
to include grid point motion. 
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In contrast to the dynamic or implicit temporal treatment of grid point 
motion, the static methods offer a great amount of simplicity, efficiency and 
spatial control at the expense of losing the accurate tracking of severe dis- 
turbances. The corrective tracking measures taken are typically to either use 
a smaller time step or preferably to require a broader band of resolution. 

With the broader band, the idea is that the disturbance will still appear in 
the high resolution region at the end of the time step. Methods which lead to 
such breadth typically come from grid smoothness forces and from curvature 
clustering on the monitor surface. 

The static methods also tend to offer more numerical stability for the 
class of problems where the broad-banded resolution provides an adequate buf- 
fer for the containment of the disturbance. In the application, either the 
grid velocities are employed with backward differences in time or the solution 
is simply interpolated onto the new grid point locations in what is known as 
remapping step. While both are commonly used, the remapping approach is more 
prevalent , particularly in cases where the steady-state convergence is a prime 
element. In the steady-state cases, the movement may start with direct inter- 
lacing between the PDE-solver and the grid generator and then progress to few- 
er and fewer applications of grid movement until movement is stopped altoge^h 
er. An example of the interlaced approach is given by Palmerio and Dervleux 
[94], In addition, they considered specified grids as in Eq. 83 and a form of 
time filtering similar to that of Winkler, Mihalas and Norman [Ml]. 

With trapezoidal finite elements in one space dimension and time, Davis 
and Flaherty [95] and Flaherty et al. [96] employed a static grid generation 
scheme for a PDE-solver that was well-adjusted for temporal evolution. The 
static grid was generated from data at the explicit time level n for use at 
the implicit time level n+1 . In a sense, this represents a zeroth-order for- 
ward extrapolation of the grid in time. Because of the extrapolation, the 
tracking possibilities are absent. In contrast, Sanz-Serna and Christie [97] 
and Blom, Sanz-Serna and Verwer [98] consider a predictor step for the appll 
cation of static grid generation. The essence of their idea is to apply the 
PDE-solver to get provisional solution values at n+1 and from those values to 
generate the grid at the implicit level n+1. Then with the grid determined at 
n+1. they get the actual solution at n+1. Although their development was one- 
dimensional in space, the idea of first predicting the implicit-level adaptive 
data is attractive for any application in any number of spatial dimensions. 
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The extra cost only amounts to an extra application of the PDE-solver. The 
benefits are the capability to accurately track rapidly moving severe dis- 
turbances by using static adaptive grid generators which are known to produce 
high quality versatile grids. 


CONCLUSION 

Adaptive grid generation is essential when rapidly varying solutions to 
partial differential equations are to be simulated in an orderly fashion. The 
fundamental character of the topic is that the necessary local resolution is 
dynamically provided while the regular ordering of points is preserved. With 
the preserved order, many of the best partial differential equation solvers 
are available and are easily formulated. The main spirit behind the various 
developments is the creation of effective grid movement strategies that can be 
coupled into a wide variety of PDE-solvers. 

In our discussion, we have attempted to develop the topic of adaptive 
grid generation in a somewhat structured and coherent manner. The emphasis 
was on the basic concepts and the interrelationship between the various con- 
sequent methods. To maintain a general perspective, basic principles were 
considered from various viewpoints. This started with a dozen basic ways to 
view equidi 3 tribution in one dimension. It recurred again when the use of 
prescribed distributions and coefficients in the weights were explored. Upon 
a foundation of concepts developed in one dimension, the various higher-di 
mensional methods were then developed. The discussion closed with a con- 
sideration of the temporal coupling of PDE- 3 olvers and grid generators. 
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A general purpose adaptive grid method will be presented that has the ability to generate time 
accurate grids for a wide variety of problems. The solution method consists of three parts, a 
grid movement scheme; a PDE solver; and a temporal coupling routine that links the dynamic 
grid and the PDE solver. The ability of the basic scheme to perform time accurate adaptive 
computations has been previously established [1|. Here, we will present a new innovation for 
grid control, address efficiency improvements in the temporal coupling and issues pertaining to 
the transfer of solution data between successive grids, and demonstrate the adaptive solution 
method on a supersonic flow problem. 

The basic grid movement scheme employs a “monitor surface” formed from the solution data 
to identify regions requiring further grid resolution. The grid points are repositioned on a 
curve by curve basis by requiring that a weight function bn equally distributed over each 
curve. Including the monitor surface gradient and normal curvature in the weight function 
ensures that both the gradients and the transition regions of the solution are resolved, ihe 
grid movement, is interwoven with a smoothing operation to ensure grid smoothness an o 
alleviate grid skewness. The scheme also contains a means of eliminating outlier values m 
the adaptive data and a grid control that can enforce a prescribed minimum grid cell size. 

A new grid control will be presented that can accurately track multiple solution features. 
At present, most adaptive grid methods can only track a single solution feature or at best a 
linear combination of multiple features. Treating the adaptive data as a aca/ar function can 
result in an inaccurate grid if the solution features should merge [2]. This limitation can be 
overcome by combining the individual features to be tracked into a monitor surface which is a 
vector function [2l. With respect to the vector monitor surface, examples of the improved grid 
resolution it provides and the required modifications for computing its curvature properties 
will be discussed. 


A simple predictor corrector scheme is used to couple the adaptive grid to the PDE solver. 
The scheme treats the time integration as a series of initial value problems m which the solution 
is first advanced on the existing grid for the purposes of obtaining a new grid, the solution data 
is then interpolated from the old grid to the new grid, and last, the solution is recomputed 
on the new grid. This approach creates a time accurate grid because the grid does not lag 
the solution in time. Furthermore, grid velocity terms are not required because the solution 
is computed on a static grid. Techniques will be presented that reduce the computational 
time consumed by the algorithm described in earlier work [1]. In addition, tests to determine 
the effect on the solution of using a conservative data transfer method, as opposed to the 
bilinear interpolation used at present, will be summarized. 

The capabilities of the adaptive method will be demonstrated by computing a time accurate 
solution for the inviscid unsteady flow field created by a shock vortex interaction. The shock 
vortex problem consists of an initially planar shock wave marching toward and eventually over 
a solid core vortex lying a short distance upstream. Here, the objective is to generate a grid 
that will track the shock wave, the vortex, and the sound wave emitted by the interaction; 
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these three features have substantially different magnitudes, and locations. Thus, the shock 
vortex problem is a particularly good test of an adaptive method because it requires a grid 
with locally high resolution whose location changes with time to properly capture both the 
severe behavior of the shock wave and the much more gentle features of the sound wave and 
the vortex. In this study we will focus on problems in which the shock wave is propagating 
at a Mach number slightly greater than 1 and thereby results in more localized behavior that 
must be refined, as compared to our earlier work in which the shock wave was propagating 
at Mach 3 [1,2]. A preliminary result in which the shock wave is propagating at Mach 1.1 is 
contained in Fig. 1. 
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R 1 For the shock vortex problem: problem definition (a), grid and pressure contour plots 
before (b) and after (c) interaction for a grid adapted to the shock wave and the vortex 
using the vector monitor surface. 




